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Abstract 


One  of  the  challenges  associated  with  defending  against  ballistic  missiles  is  to 
discriminate  the  object  of  interest  among  multiple  closely-spaced  objects  (CSOs)  that  travel  on  a 
ballistic  trajectory.  The  discrimination  process  typically  involves  the  identification  and  tracking 
of  the  object  of  interest.  One  approach  that  can  improve  discrimination  perfonnance  is  to 
employ  multiple  sensors.  Multiple-sensor  correlation  and  discrimination  involves  the  integration 
of  sensor  measurements  collected  from  terrestrial  and  on-orbit  sensors  to  improve  the  likelihood 
of  identifying  and  tracking  an  object  of  interest  within  the  CSOs.  This  report  describes  the 
development  of  the  algorithms  necessary  for  fusing  sensor  measurement  data  obtained  from 
multiple,  dissimilar  sensors  in  order  to  improve  the  likelihood  of  identifying  and  tracking  an 
object  of  interest  within  closely-spaced  objects  traveling  on  a  ballistic  trajectory.  The  algorithms 
utilize  a  target  object  map  (TOM)  that  is  created  using  multiple  sensor  measurements  for 
correlation.  The  object  of  interest  is  then  selected  using  a  probability-based  Dempster-Shafer 
discrimination  algorithm.  To  examine  the  performance  of  these  algorithms,  a  simulation 
environment  was  developed.  It  included  relevant  characteristics  of  the  sensors  in  the 
discrimination  system,  a  modeling  process  for  the  ballistic  trajectories  of  the  CSOs,  and  a 
decision-making  module  containing  the  algorithms  for  handling  the  sensor  returns  and 
correlating  and  discriminating  the  object  of  interest.  This  simulation  environment  was  iterated  to 
assess  the  effect  of  varying  radar  angular  and  range  resolution  and  discrimination  characteristic 
distribution  overlap  on  discrimination  performance.  This  project  found  there  was  not  a  strong 
correlation  between  changing  radar  resolution  by  +/- 1 0%  from  the  nominal  values  and 
probability  of  discrimination.  However,  it  was  noted  that  a  35  percent  increase  in  track 
uncertainty  for  an  increase  of  0.1  degrees  in  angular  resolution  and  a  six  percent  increase  in  zero 
effort  miss  distance  for  the  same  change  in  angular  resolution.  This  project  also  noted  a  strong 
inverse  quadratic  relationship  between  the  characteristic  overlap,  which  is  determined  by  the 
amount  of  overlap  in  the  probability  of  observation  for  different  classes  of  objects,  and 
discrimination  performance. 


Keywords:  (1)  closely-spaced  ballistic  objects,  (2)  multiple,  dissimilar  sensor  correlation,  (3) 
target  object  map,  (4)  Dempster-Shafer  discrimination  algorithms,  (5)  simulation  environment 
for  various  operational  scenarios,  (6)  sensitivity  analysis 
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Preface 


Ballistic  missiles  have  long  been  a  challenging  tactical,  operational,  and  strategic  asset  as 
well  as  problem  for  militaries  around  the  world.  The  strategic  stability  of  the  Cold  War  was 
largely  predicated  on  the  invulnerability  of  nuclear-tipped  land  and  submarine  launched 
intercontinental  ballistic  missiles.1  However,  as  the  technologies  for  these  systems  have 
matriculated  from  major  strategic  powers  to  smaller,  less  predictable  rogue  states,  a  new  tactical 
and  strategic  threat  has  emerged."  The  performance  of  the  Patriot  system  in  the  Gulf  War, 
shooting  down  51  of  88  potential  threats  at  a  cost  of  157  expended  interceptors3,  highlighted  the 
need  for  better  missile  defense  systems  to  protect  theater  forces  from  ballistic  missile  threats, 
while  the  possible  threat  of  a  rogue  missile  attack  on  the  United  States  homeland  motivates  a 
system  capable  of  intercepting  intercontinental  threats4.  These  systems  are  necessary  not  merely 
to  defend  America  and  its  allies  from  potential  attack  but  render  ineffective  the  threat  of  potential 
attack.  In  this  manner,  “ballistic  missile  defense  is  not  simply  a  shield  but  an  enabler  of  U.S. 
action.”5  To  fail  to  develop  flexible,  effective  missile  defense  systems  would  be  to  surrender  the 
strategic  imperative  around  the  world,  a  fate  that  should  be  avoided  by  the  United  States. 
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1  Introduction 


Discrimination,  at  its  heart,  is  nothing  more  than  decision  making.  It  is  defined  as  the 
ability  to  make  fine  distinctions.  Discrimination  in  a  technical  context  is  then  the  synthesis  of  the 
given  information  in  a  coherent  manner  in  order  to  “make  a  fine  distinction.”  Various  factors 
convolute  and  confuse  this  process.  For  example,  the  information  utilized  for  making  this 
decision  is  almost  always  imperfect  because  of  sensor  error,  and  in  many  cases  information  from 
different  sources  can  support  different  decisions.  In  some  scenarios  such  as  military  applications, 
the  situation  often  demands  dynamic  and  responsive  decisions  due  to  challenging  environments. 
In  these  situations  a  wrong  decision  may  lead  to  dire  consequences. 

While  discrimination  in  general  has  a  variety  of  applications,  this  project  is  focused  on 
discrimination  in  the  context  of  Ballistic  Missile  Defense  (BMD).  In  a  typical  BMD  system,  a  set 
of  sensors  is  responsible  for  observing  a  complex  of  objects,  selecting  a  target  to  engage,  and 
intercepting  that  target.6  A  multiple-sensor  discrimination  system  then  utilizes  more  than  one 
sensor  to  accomplish  the  discrimination  task.  In  this  context,  the  target  object  the  system  needs  to 
select  is  the  Re-entry  Vehicle  (RV),  which  is  the  lethal  object  in  the  complex.  A  ballistic  missile 
launch  will  also  yield  a  spent  booster  section,  known  as  the  tank,  natural  debris  from  the 
deployment  of  objects,  and  countenneasures  to  deceive  a  defense  system.7  Multiple  sensors  offer 

the  ability  to  observe  different  sets  of  phenomena  and,  when  properly  synthesized,  can  improve 

8 

discrimination  perfonnance. 

This  scenario  is  complicated  by  three  principal  factors.  First,  the  physics  of 
exoatmospheric  ballistic  trajectories  present  difficult  viewing  conditions  for  each  sensor  that 
makes  sensor  information  difficult  to  correlate.  On-orbit  objects  velocities  of  approximately 
seven  km/s  can  significantly  limit  the  observation  window  of  the  sensors  on  the  ground  while 
closing  velocities  in  excess  of  10  km/s  can  limit  the  observation  window  for  interceptors.9 
Second,  Closely-Spaced  Objects  (CSO)  present  significant  difficulties  for  sensor  observation  and 
multiple-sensor  correlation  as  sensor  measurement  uncertainty  makes  high  fidelity  position 
determination  difficult.  “Closely-spaced,”  in  this  context,  can  be  regarded  as  objects  that  are 
spaced  closer  than  the  measurement  uncertainty  band  for  one  or  multiple  sensors  within  the 
observation  architecture.  Resolving  targets  from  one  sensor  image  to  another  can  prove  difficult 
under  these  conditions.  Finally,  the  sensors,  such  as  radar  and  optical  sensors,  are  corrupted  by 
noise  and  bias  making  discrimination  still  more  difficult. 

1.1  Problem  Statement 

This  report  describes  the  development  of  the  necessary  algorithms  for  the  multiple  sensor 
correlation  and  discrimination  of  CSOs 10,11  and  the  evaluation  of  their  perfonnance  in  a  relevant 
simulation  environment.  These  algorithms  handle  the  sensor  returns,  process  them  into  useful 
metrics,  and  assess  the  probability  that  a  given  object  is  the  object  of  interest  from  these  metrics. 
To  examine  the  performance  of  these  algorithms,  a  simulation  environment  was  developed.  It 
included  relevant  characteristics  of  the  sensors  in  the  discrimination  system,  a  modeling  process 
for  the  ballistic  trajectories  of  the  CSOs  in  the  target  complex,  and  a  decision-making  module 
containing  the  algorithms  for  handling  the  sensor  returns  and  correlating  and  discriminating  the 
object  of  interest. 
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The  measure  of  discrimination  performance  on  a  ballistic  missile  target  complex  for  a 
multi-sensor  system  is  the  probability  of  the  ballistic  missile  defense  system  properly 
discriminating  the  RV  from  other  objects  in  the  target  complex.  This  probability  is  known  as  the 
probability  of  discrimination,  P A  and  depends  upon  the  system’s  ability  to  accurately  track  and 
correlate  each  object  in  addition  to  correctly  discriminating  each  object  and  selecting  the  RV. 
The  Figure  of  Merit  (FOM)  to  assess  the  system’s  correlation  perfonnance  is  the  probability  of 
correlation,  Pc. 

1.2  Nominal  System  Architecture 

A  multiple-sensor  discrimination  system  will  be  comprised  of  at  least  two  sensors.  As 
dictated  by  the  need  to  conduct  early  warning,  search,  sensor  cueing,  tracking,  discrimination, 
and  homing,  a  system  can  have  additional  sensors  as  necessary  to  meet  these  needs.  However, 
this  research  is  primarily  concerned  with  analyzing  the  discrimination  performance  of  the 
system.  To  this  end,  only  the  sensor  components  essential  to  discrimination  will  be  considered. 
The  nominal  system  consisted  of  a  sea-based  radar  system  and  an  optical  sensor  deployed  in  a 
Kill  Vehicle  (KV)  to  intercept  the  target.  This  architecture  is  depicted  in  Figure  1. 


Figure  1 .  Nominal  system  architecture 

The  radar  system  was  responsible  for  tracking  the  complex,  providing  data  to  calculate  a  launch 
solution  to  intercept  the  target  complex,  and  providing  discrimination  information.  The  KV 
possessed  an  onboard  optical  sensor  for  target  acquisition,  tracking,  and  discrimination.  The 
system  fused  the  tracks  and  discrimination  information  from  the  two  sensors  to  select  a  target. 
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2  Simulation  Environment  Development 


Figure  2  depicts  the  steps  performed  in  the  simulation  environment  and  the  requisite 
elements  for  modeling  that  particular  step.  The  simulation  depicted  the  launch  of  a  threat, 
followed  by  a  radar  tracking  the  threat,  the  launch  of  an  interceptor,  and  finally  the  tenninal 
acquisition  by  the  Kill  Vehicle  (KV). 
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Figure  2.  Simulation  environment  developed  to  support  a  concept  of  operation 

2,1  Target  Complex  Trajectory  Modeling 

Two-body  orbit  dynamics  was  utilized  to  generate  ballistic  trajectories  between  the 
desired  burnout  and  reentry  points.  The  dispersion  of  object  trajectories  within  the  target 
complex  was  achieved  by  varying  the  velocity  of  each  object  at  the  burnout  points.  The 
generated  trajectories  served  as  the  "truth"  for  the  sensor  modeling  process.  The  System  Tool  Kit 
(STK)  was  utilized  to  calculate  the  trajectories  of  these  objects.  This  interface  is  depicted  in 
Figure  3.  The  interface  allowed  STK  to  be  controlled  iteratively  and  autonomously  by 
MATLAB.  Functions  such  as  simulation  control,  tracking,  correlation,  and  discrimination  took 
place  within  MATLAB.  The  necessary  information  and  commands  were  then  passed  to  STK  to 
generate  truth  and  estimated  trajectories  and  calculate  the  interceptor  flight  path.  The  entirety  of 
the  interface  code  is  contained  in  Appendix  IV.  The  complete  interceptor  trajectory,  from  initial 
launch  to  target  intercept,  was  considered  in  the  concept  of  operation.  However,  only  the 
exoatmospheric  trajectory  of  the  target  complex  was  modeled,  starting  immediately  after 
complex  deployment  until  re-entry. 
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MATLAB: 

Simulation  Control 
Radar  Tracking 
Correlation 
Discrimination 


MATLAB-STK  Interface 


Figure  3.  MATLAB  -  STK  interface 

2.2  Radar  Observation  Modeling 

A  radar  system  takes  measurements  of  azimuth,  elevation,  and  range  for  each  object  it 
detects,  and  calculates  the  associated  rates  of  these  measurements.  To  model  a  radar  observation, 
measurement  error  was  added  to  the  truth  measurements  of  the  target  complex.  Measurement 
error  consisted  of  two  sources,  accuracy  and  precision.  Accuracy  was  defined  as  the  degree  to 
which  the  sensor  was  able  to  measure  truth,  while  precision  was  the  uncertainty  associated  with  a 
measurement.  Accuracy  was  reflected  by  the  bias  of  the  radar  system.  It  was  assumed  that  the 
bias  was  constant  for  the  radar  system,  meaning  the  same  error  in  azimuth,  elevation,  and  range 
was  added  to  each  detected  target.  Precision  was  reflected  by  the  random  noise  in  each 
subsequent  measurement.  The  noise  was  assumed  to  be  Gaussian  with  a  mean  of  zero.  The 
standard  deviation  of  the  noise  distribution  was  based  upon  the  radar  range  and  angular 
resolution.  The  measurements  were  then  given  by  Equation  (  1  ),  where  z*  is  the  noisy 
measurement,  jc*  is  the  truth  measurement,  w*  is  the  noise  at  a  discrete  moment  in  time,  and  bk  is 
the  bias. 


Zk  =  xk  +  wk+bk.  (1) 

The  radar  system’s  ultimate  task  was  to  track  the  target  complex  in  order  to  estimate  its 
ballistic  trajectory.  Due  to  the  dynamic  nature  of  the  trajectory,  a  Kalman  Filter  (KF)  was 
implemented  to  estimate  the  position  and  velocity  of  the  target  in  the  topocentric-horizon 
coordinate  system,  also  known  as  the  South,  East,  Zenith  (SEZ)  coordinate  frame.  The 
MATLAB  code  for  the  KF  is  contained  in  Appendix  III.  This  coordinate  frame  was  centered  at 
the  radar  site.  Track  information  was  then  converted  into  the  Earth-Centered  Inertial  (ECI) 
coordinate  system  for  orbit  ephemeris  determination.  The  SEZ  and  ECI  coordinate  frames  are 
depicted  in  Figure  4. 
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Figure  4.  SEZ  and  ECI  coordinate  frames 

2.3  Optical  Observation  Modeling 

Unlike  radar,  an  optical  sensor  is  only  capable  of  capturing  two-dimensional  position 
information  within  the  focal  plane's  field  of  view.  In  order  to  model  an  optical  observation,  the 
truth  complex  had  to  be  projected  onto  the  field  of  view  of  the  KV.  First,  the  position  of  the 
target  complex  was  described  in  the  KV-centered  coordinate  system,  depicted  in  Figure  5,  where 
the  p-axis  is  along  the  boresight  from  the  KV  to  the  RV  and  the  m-  and  n-axes  are  the 
perpendicular  axes. 
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This  description  required  knowledge  of  the  truth  position  of  the  KV  and  the  target 
complex  at  each  measurement  instance.  The  coordinate  transformation  from  ECI  to  p-m-n  was 
then  accomplished  using  Euler  angles.  To  transform  from  the  ECI  coordinates  into  the  p-m-n 
frame,  the  relative  position  of  the  target  complex  from  the  KV  was  calculated  in  ECI,  using 
Equation  (  2  ).  P  is  the  relative  position  of  the  target  complex  from  the  KV,  R,  is  the  absolute 
position  of  the  target  complex,  and  R,  is  the  absolute  position  of  the  KV.  A  similar  process  was 
utilized  to  calculate  the  relative  velocity  of  the  target  complex. 

-  =  !k~!k  (2) 

With  the  relative  position  and  velocity  known,  the  ECI  coordinates  could  then  be 
transformed  into  p-m-n  using  a  Direction  Cosine  Matrix  (DCM),  which  represents  rotation  about 
the  3  (K)  axis  followed  by  rotation  about  the  2  (J)  axis.  First,  the  ECI  frame  was  rotated  about 
the  K-axis  by  a,  given  by  Equation  (  3  ),  where  a  and  b  are  the  I  and  J  components  of  the 
relative  position,  respectively. 


Next,  the  coordinate  frame  was  rotated  around  the  J-axis  by  P,  given  by  Equation  (  4  ), 
where  c  is  the  K  component  of  the  relative  position,  and  d  is  the  magnitude  of  the  I-J 
components. 

p  =  an-1  (-j)  (4) 

The  DCM  was  then  assembled  from  the  constituent  components,  as  shown  in  Equations 
(  5  )  -  (  7  ),  where  D  is  the  DCM  used  to  rotate  from  ECI  to  p-m-n. 
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With  the  DCM  known,  the  p-m-n  position  could  be  calculated  via  Equation  (  8  ),  where 
P’  is  the  position  of  the  target  complex  in  p-m-n  coordinates.  Velocity  could  be  found  in  a 
similar  manner. 


With  three-dimensional  position  and  velocity  known,  range,  azimuth,  and  elevation  were 
calculated  based  upon  the  Cartesian  p-m-n  elements.  The  optical  observations  were  assumed  to 
be  the  truth  observations. 
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3  Correlation  and  Discrimination  Algorithms 


3.1  Target  Object  Map  Generation  and  Correlation 

With  the  state  estimates  determined  by  the  radar  and  subsequently  observed  by  the  KV 
optical  sensor,  the  infonnation  from  these  two  sources  needed  to  be  correlated.  To  accomplish 
this,  the  system  must  be  able  to  associate  the  different  sensor  returns  from  the  same  object.  A 
method  known  as  anchor  Target  Object  Map  (TOM)  is  used.  A  TOM  is  essentially  a  picture  of 
the  target  complex  at  a  discrete  time  generated  by  a  sensor.  The  TOM  is  a  matrix  containing  the 
coordinates  of  the  location  of  the  center  of  mass  of  each  of  the  detected  objects  in  the  target 
complex  in  a  reference  coordinate  frame.  TOMs  from  different  sensors  can  then  be  compared 
provided  that  each  TOM  reflects  similar  viewing  geometry  and  the  same  instance  in  time. 
Because  the  sensors  within  the  architecture  rarely  observe  a  target  at  the  same  instance  from  the 
same  location,  a  method  was  developed  to  generate  similar  TOMs  both  temporally  and 
geometrically.  The  MATLAB  code  utilized  for  this  method  is  contained  in  Appendix  V. 

The  radar  TOM  was  generated  using  the  tracks  from  the  radar  observation.  First,  the 
target  complex  track  was  projected  from  the  time  of  minimum  track  uncertainty,  defined  as  the 
nonnal  of  the  variance  elements  of  the  track  in  the  KF  estimates,  to  the  desired  time.  With  the 
anticipated  position  known,  the  coordinates  of  the  target  complex  were  then  rotated  and 
projected  into  the  focal  plane  of  the  KV,  similar  to  the  procedure  described  in  Section  2.3.  The 
optical  TOM  could  be  obtained  directly  from  the  KV  observations  since  the  radar  TOM  is 
projected  to  the  orientation  and  time  during  which  the  KV  is  observing  the  complex. 

With  a  TOM  generated  for  each  sensor,  the  TOMs  had  to  be  correlated  as  depicted  in 
Figure  6.  First,  the  sensor  bias  was  eliminated  by  “anchoring”  each  TOM  to  a  reference  object. 
The  reference  object  is  one  readily  identifiable  by  both  sensors.  The  anchor  was  selected  using 
the  probabilistic  discrimination  algorithm  described  in  Section  3.2.  Each  sensor  independently 
performed  discrimination  to  select  its  own  anchor  for  its  TOM.  After  selecting  an  anchor,  the 
anchor  was  set  as  the  origin  of  the  TOM  by  subtracting  its  coordinates  from  the  coordinates  of  all 
the  objects. 
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RadarTOM 


Optical  Sensor 


COMBINED  TOM 


Figure  6.  TOM  anchoring  and  correlation  processes 

The  two  TOMs  were  then  superimposed.  Due  to  sensor  uncertainty,  the  observed  location 
of  each  object  may  differ  between  the  two  sensors.  As  the  spacing  between  the  objects  decreases, 
the  relative  arrangement  of  the  target  complex  may  vary  between  two  TOMs.  The  differing  in 
target  position  makes  the  association  of  objects  difficult.  To  correctly  associate  the  targets,  each 
possible  configuration  of  associations  was  tested.  An  example  of  this  correlation  process  for  four 
active  tracks  from  the  radar  and  the  optical  sensor  is  depicted  in  Figure  7.  Pairings  between 
objects  are  depicted  by  the  green  lines  between  objects. 
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Figure  7.  Possible  track  correlation  pairings  for  four  objects 

For  each  arrangement,  a  correlation  coefficient,  r,  was  calculated  using  Equation  (  9  ),  where  A 
and  B  are  the  matrices  containing  the  coordinates  within  the  TOM  and  m  and  n  are  the  number  of 
observations  within  A  and  B ,  respectively.  Maximizing  r  corresponded  to  minimizing  the 
Euclidean  distance  between  the  pairings  in  the  TOM.  A  nominal  TOM  A  matrix  is  depicted  in 
Table  1.  Note  that  the  track  numbers  (shaded  in  red)  are  only  for  reference  and  are  not  included 
for  the  purposes  of  calculating  r.  In  this  case,  track  2  was  identified  as  the  anchor  and  as  such  is 
set  as  the  origin  (0,  0).  A  and  B  are  the  mean  values  of  the  A  and  B  matrices. 

Table  1.  Example  TOMA  matrix 


Track 

Azimuth  (rad) 

Elevation  (rad) 

1 

-0.014 

0.015 

2 

0 

0 

3 

0.030 

-0.045 

4 

0.039 

-0.049 
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The  track  correlation  pairing  generating  the  greatest  r  was  then  selected  and  the  respective  tracks 
are  associated. 

3.2  Discrimination  Algorithm 

The  discrimination  algorithm  was  responsible  for  selecting  a  target  from  the  sensors’ 
discrimination  information.  The  algorithm  used  Dempster-Shafer  decision  logic,  which  is 
explained  in  Section  3.2.1  -  Section  3.2.4.  The  MATLAB  implementation  of  this  algorithm  is 
contained  in  Appendix  VI. 

3.2.1  Dempster-Shafer  Decision  Logic 

As  previously  stated,  discrimination  is  the  process  of  transforming  measurements  into  a 
decision.  In  the  context  of  this  project,  the  discrimination  algorithm  is  utilized  to  select  the  RV 
given  a  set  of  measurements.  Dempster-Shafer  logic  was  employed  within  this  algorithm,  as  the 
Dempster-Shafer  algorithm  “is  useful  when  the  sensors  contributing  information  cannot  associate 
a  100  percent  probability  of  certainty  to  their  output  decisions.”14  The  Dempster-Shafer  method 
captures  and  combines  the  certainty  each  sensor  possess  in  its  object  discrimination 
capabilities.15  Uncertainty  in  sensor  measurement  is  of  particular  concern  when  discriminating 
targets  on  a  ballistic  trajectory.  The  extreme  geometry  of  the  problem,  the  associated  temporal 
constraints,  and  ambiguous  and  imperfect  sources  of  information  all  contribute  to  the  uncertainty 
in  selecting  the  RV  from  the  target  complex. 

Dempster-Shafer  involves  the  combination  of  data  from  a  set  of  sensors.  Each  sensor 
observes  a  certain  set  of  characteristics  and  generates  returns  from  these  observations.  From 
these  observations,  the  probability  of  each  measurement  resulting  from  a  certain  type  of  object, 
or  object  class,  is  calculated.  This  probability  of  each  class  is  used  to  calculate  a  belief  in  a  given 
hypothesis,  known  as  an  Evidence  Measure  (EM)  that  reflects  the  relative  certainty  of  selecting  a 
certain  hypothesis.  The  hypotheses  are  broken  into  two  categories:  singleton  and  combination. 
Singleton  hypotheses  categorize  the  observation  into  one  object,  such  as  an  RV.  Combination 
hypotheses  categorize  the  observation  into  multiple  objects,  e.g.  RV  or  decoy.  An  EM  of  unity 
indicates  absolute  certainty  in  a  certain  hypothesis.  The  EMs  from  every  sensor  are  then 
combined  utilizing  Dempster’s  rule  of  combination.  The  hypothesis  with  the  greatest  amount  of 
accumulated  EM  is  then  selected. 

3.2.2  Calculating  Evidence  Measures 

In  order  to  calculate  the  EM  for  each  observation,  the  distribution  of  measurements  must 
be  known.  This  project  assumes  this  distribution  to  be  Gaussian,  with  the  expected  value  for  a 
given  class  to  be  the  mean  of  the  distribution,  p.16  The  sensor  measurement  uncertainty 
determines  the  standard  deviation  of  the  distribution,  a.  From  the  Gaussian  distribution,  the 
probability  of  a  given  measurement  stemming  from  a  given  object  can  be  calculated.  Figure  8 
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depicts  an  observation  and  the  calculation  of  the  probability  of  measurement  for  a  three-class 
hypothesis  example.  This  observed  characteristic  can  be  anything  the  sensor  is  capable  of 
measuring,  such  as  the  size  or  a  dynamic  characteristic  of  the  object.  This  project  assumed  the 
sensors  are  equally  capable  of  measuring  any  class  of  object.  As  such,  the  characteristic 
distributions  are  symmetric  but  possess  different  expected  values.  Table  2  shows  the  calculated 
probabilities  from  these  distributions. 


Figure  8.  Example  calculation  of  measurement  probabilities  for  three-class  hypotheses 
Table  2.  Example  calculated  probabilities 


P(XRV) 

0.629 

P(X  Tank) 

0.021 

P(XOther) 

~0 

A  given  set  of  EMs  must  then  be  assigned  to  support  a  series  of  hypotheses.  The  set  of  all 
the  possible  hypotheses  is  referred  to  as  the  frame  of  discernment,  denoted  by  0.  Any  EM  can  be 
assigned  to  any  of  the  propositions,  aj,  a2,  a„,  or  any  of  the  possible  unions,  e.g.  ai  U  a2,  or 

finally  to  0  itself.  A  given  EM  can  only  be  assigned  to  the  extent  of  the  confidence  in  the 
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measure.  For  example,  if  a  sensor  is  80  percent  confident  a  measurement  indicates  ai,  an  EM  of 
0.8  is  assigned  to  a j.  The  remainder,  0.2,  is  assigned  to  0,  which  represents  all  the  possible 
propositions  simultaneously,  or  uncertainty. 

3.2.3  Evidence  Measure  Combination 

With  the  proper  EM  assigned  to  each  hypothesis  for  each  sensor  observation,  the  EMs 
can  be  combined  to  generate  the  total  support  for  each  hypothesis  using  Dempster’s  rule. 
Dempster’s  rule  is  implemented  by  fonning  a  matrix  with  the  elements  to  be  combined  in  the 
first  column  and  first  row.  The  combinations  to  be  formed  are  those  hypotheses  which  have  an 
intersection  that  exists.  For  example,  hypotheses  ai  U  as  and  ai  U  a4  intersection  is  03. 

Combining  the  EMs  for  these  two  hypotheses  would  yield  evidence  for  a 3. 17  The  matrix  elements 
and  subsequent  combinations  are  made  by  multiplying  the  column  element  in  the  first  column 
with  the  row  element  in  the  first  row.  An  example  for  two  sensors,  sensor  A  and  sensor  B,  is 
shown  in  Table  3.  For  additional  sensors,  the  above  method  is  repeated  for  each  subsequent 
sensor. 


1 8 

Table  3.  Application  of  Dempster's  Rule 


mB(a3  U  a 4)  =0.7 

mB(0)  =  0.3 

mA(0)  =  0.4 

m(a3  U  a 4)  =  0.28 

3 

II 

* — 4 

mA(aI  U  a3)  =0.6 

m(a3)  =  0.42 

m(a3  U  a3)  =  0.18 

Sensor  A  Evidence  Measure 
Sensor  B  Evidence  Measure 
Combined  Evidence  Measure 


A  special  case  of  combination  occurs  when  measures  support  mutually  exclusive  hypotheses,  for 
example  ai  and  ct2  U  as.  In  this  case,  the  product  of  these  two  EMs  is  assigned  to  the  empty  set, 

<f>,  which  represents  measures  that  support  different  and  exclusive  hypotheses.  (f>,  also  known  as 
conflict,  is  the  complement  of  0,  which  represents  a  measure  that  could  represent  any  and  all 
hypotheses.  The  empty  set  (j>  is  only  pertinent  in  the  combination  of  EMs,  as  a  single  measure  can 
be  uncertain  in  what  hypothesis  it  supports  but  cannot  conflict  itself. 

3.2.4  Resolving  Conflict  and  Uncertainty 

While  knowledge  of  a  degree  of  conflict  and  uncertainty  in  some  applications  may  prove 
useful,  they  are  not  consistent  with  BMD  discrimination,  hi  this  context,  uncertainty  and  conflict 
must  be  eliminated  and  assigned  to  support  a  hypothesis.  To  do  so,  one  must  calculate  the 
“pignistic  probability”  of  a  certain  hypothesis  from  the  given  EMs.  Pignistic  is  derived  from  the 
Latin  term  pignus,  meaning  a  bet,  and  is  indicative  of  the  nature  of  Pignistic  Probabilities  (PP).19 
PPs  are  not  analytical  probabilities  in  the  sense  that  a  PP  of  0.5  would  not  indicate  an  event  is 
expected  to  occur  50  times  in  100  trials.  Rather,  they  are  meant  to  guide  a  “bet,”  or  decision,  to  a 
particular  hypothesis.  Their  absolute  value  is  irrelevant,  only  their  relative  magnitudes  matter. 

For  example,  consider  the  two  PPs  for  hypotheses  a  and  b ,  where  PPa=  0.75  and  PPb  =  0.25.  The 
pignistic  probability  of  a  does  not  indicate  three  times  greater  confidence  in  a,  merely  that  a 
should  be  selected  over  b. 

To  calculate  the  PP  for  a  given  hypothesis,  conflict  is  first  addressed.  To  eliminate 
conflict,  the  EMs  for  the  non-empty  set  hypotheses  are  normalized  by  a  factor,  K,  given  by 
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Equation  (  10  ).  K  is  a  measure  of  the  total  amount  of  conflict  between  two  measurements, 
where  a  if  of  unity  indicates  no  conflict  and  an  undefined  K  indicates  complete  conflict  between 
the  measurements,  causing  the  sum  given  by  Dempster’s  rule  to  be  undefined. 

K-1  =  1  -  ^  [ma(ai)mi,(fc/)]  (  10  ) 

at  n  b j=0 

After  normalizing,  the  conflict  terms  are  set  to  zero.  With  conflict  eliminated,  the  EMs  for 
combination  hypotheses  and  uncertainty  are  then  equally  divided  between  their  constituent 
components.  For  example,  if  m(aiUa2)  =  0.5,  m(ai)  and  m(a2)  would  both  be  assigned  0.25.  The 
uncertainty  EM  would  then  be  divided  equally  among  all  the  hypotheses.  At  this  point,  the 
remaining  hypotheses  are  the  singleton  elements  ( aj ,  a?,  ...).  The  evidence  measure  now 
assigned  to  these  elements  is  the  PP  for  each  hypothesis.  With  conflict,  uncertainty,  and  the 
combination  hypotheses  eliminated,  the  hypothesis  with  the  highest  total  PP  is  then  selected. 
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4  Nominal  Scenario  Simulation  and  Results 


In  order  to  support  future  decisions  regarding  investment  in  and  improvement  of  the 
ballistic  missile  defense  system,  the  effect  of  varying  system  parameters  on  Pd  was  assessed 
using  a  simulation  environment  developed  as  a  part  of  this  project.  The  simulation  environment 
depicted  a  launch  event  of  a  threat,  tracking  by  a  radar  system,  and  finally  engagement  by  an 
interceptor.  The  simulation  then  assessed  the  success  or  failure  of  the  engagement.  The  nominal 
scenario  and  one  iteration  of  the  simulation  are  described  in  Section  4.1  through  Section  4.4. 

4.1  Nominal  Simulation  Scenario 

To  assess  the  created  simulation  environment,  a  nominal  scenario  was  developed  that  was 
indicative  of  a  possible  engagement  of  a  closely-spaced  target  complex.  The  scenario  involved 
the  launch  of  a  ballistic  missile  from  the  Reagan  Test  Site  in  the  Kwajalein  Atoll  with  a 
destination  of  Seattle,  WA.  The  missile  was  assumed  to  deploy  the  target  complex  at  an  altitude 
of  100  km  directly  above  the  Reagan  Test  Site  and  re-enter  the  atmosphere  at  an  altitude  of  100 
km  directly  above  Seattle.  The  target  complex  consisted  of  a  RV,  a  tank,  and  two  decoys.  The 
trajectories  of  the  complex  objects  are  shown  in  Figure  9,  with  the  RV,  tank,  and  decoys 
trajectories  shown  in  red,  orange,  and  yellow,  respectively. 
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Figure  9.  Nominal  scenario  ballistic  trajectories 

The  defensive  system  consisted  of  a  sea-based  X-band  radar  (XI)  north  of  Hawaii  and  an 
interceptor  launched  from  Vandenberg  AFB,  CA,  that  deployed  the  KV.  The  interceptor  was 
assumed  to  have  a  range  of  2500  km,  while  XI  was  assumed  to  have  a  range  of  4500  km.  The 
interceptor  attempted  to  engage  the  target  at  an  altitude  of  250  km.  The  KV  was  assumed  to 
acquire  the  target  at  a  range  of  500  km  and  tracked  the  complex  as  it  closed.  The  radar  and 
interceptor  coverage  are  depicted  in  Figure  10,  with  XI  coverage  shown  in  light  blue  and  the 
interceptor  shown  in  light  green.  In  addition,  the  portion  of  the  target  complex  trajectory  that  can 


be  tracked  by  XI  is  highlighted  in  light  blue  and  the  portion  that  is  in  interceptor  range 
highlighted  in  light  green. 
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Figure  10.  XI  and  interceptor  coverage 

Using  the  modeling  techniques  described  in  Section  1.2,  an  engagement  of  the  target  complex 
was  simulated.  The  defensive  system  was  responsible  for  tracking  the  complex  with  both  the 
radar  (XI)  and  the  KV  optical  sensor,  correlating  the  tracks,  fusing  the  discrimination 
information,  and  ultimately  selecting  a  target. 

The  following  describes  the  results  of  the  simulation  using  the  nominal  scenario 
described  above. 

4.2  Radar  Kalman  Filter  and  Tracking  Performance 

As  mentioned  in  Section  2.2,  a  Kalman  filter  was  employed  in  order  to  estimate  the  target 
complex  states  from  the  noisy  radar  measurements.  The  performance  of  the  filter  was  measured 
by  two  metrics,  the  residuals  of  the  estimates  compared  to  the  truth  measurements  and  the 
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variance  of  the  state  estimates.  Figure  1 1  illustrates  the  position  and  velocity  residuals  of  the 
estimates  for  the  RV  track  versus  the  time  after  launch. 


Time  (s) 


Figure  1 1 .  Residuals  of  position  and  velocity  estimates  for  the  RV  track 

After  an  initialization  period  (not  shown),  the  position  residual  gradually  decreased  to 
within  a  magnitude  approximately  0.5  km.  The  velocity  residuals  fall  to  within  0.03  km/s.  Both 
the  position  and  velocity  residuals  improve  as  the  target  nears  700  seconds  after  launch.  This  is  a 
result  of  the  range  to  the  target  decreasing  which  decreases  the  covariance  of  the  measurements. 
The  residuals  then  diverge  after  the  target  passes  through  the  closest  point  of  approach  as  the 
range  begins  to  increase.  This  is  of  particular  concern  for  correlation  purposes  as  larger  residuals 
and  the  associated  larger  errors  make  associating  sensor  tracks  more  difficult  since  the 
anticipated  position  of  an  object  and  its  actual  position  will  vary  by  a  larger  amount.  After 
approximately  700  seconds  after  launch,  the  range  to  the  target  begins  to  increase.  This  drives  an 
increase  in  the  residuals. 

Figure  12  depicts  the  variance  of  the  position  estimates  for  track  3.  The  variance  of  the 
estimates  exhibited  the  same  behavior  as  the  residuals,  decreasing  in  magnitude  as  the  target 
range  decreased  and  growing  as  the  range  to  the  target  increased. 
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Figure  12.  Variance  of  SEZ  position  estimates 

4.3  TOM  Correlation 

Since  the  optical  measurements  were  assumed  to  be  truth,  tracking  and  filtering  was 
unnecessary.  The  KV  TOM  was  obtained  directly  at  the  correlation  time  from  the  optical 
observations.  The  estimated  radar  states  at  this  time  were  then  rotated  and  projected  onto  the 
focal  plane  of  the  KV.  First,  an  anchor  was  selected  by  both  sensors  based  upon  their 
discrimination  information.  The  anchor  to  be  identified  was  the  tank,  as  it  is  typically  the  largest 
object  in  the  target  complex  and  most  easily  identified  by  both  sensors.  In  this  example,  both 
sensors  properly  selected  the  tank  as  the  anchor.  This  allowed  the  TOMs  to  be  properly 
anchored.  The  tracks  then  were  correlated.  Figure  7  in  Section  3.1  depicts  the  possible  pairings. 
Configuration  one  yielded  an  r  value  of  0.997,  while  configuration  two,  three,  and  four  yielded 
0.369,  -0.23,  and  0.269,  respectively.  As  such,  configuration  one  was  selected.  This  turned  out  to 
be  the  correct  correlation. 

4.4  Discrimination 

With  the  tracks  from  the  two  sensors  correlated,  discrimination  was  then  perfonned.  Each 
sensor  observed  one  discriminable  characteristic.  The  probability  distributions  for  the  radar  and 
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optical  characteristics  are  shown  in  Figure  13  and  Figure  14,  respectively.  The  shape  of  these 
distributions  is  based  upon  the  expected  value  for  the  given  characteristic  for  each  class  and  the 
standard  deviation  of  the  measurements  from  the  sensors. 


Figure  13.  Radar  characteristic  probability  distributions  for  three  classes 


Figure  14.  Optical  characteristic  probability  distributions  for  three  classes 


The  radar  and  optical  discrimination  characteristic  observations  are  contained  in  Table  4. 
Table  4.  Example  scenario  radar  and  optical  observed  discrimination  characteristics 
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Track 

Radar  Characteristic 

Optical  Characteristic 

1 

1.308 

4.094 

2 

3.234 

5.748 

3 

1.634 

5.357 

4 

2.546 

4.625 

From  these  observed  characteristics,  the  probability  of  observation  for  each  class  was  calculated 
for  each  set  of  measurements  and  is  contained  in  Table  5  and  Table  6,  respectively. 

Table  5.  Example  scenario  probability  of  observation  from  radar  characteristic  observation 


Track 

Probability  of  RV 

Probability  of  Tank 

Probability  of  Decoy 

1 

0.660 

0.003 

0.306 

2 

3.69e-5 

0.715 

0.038 

3 

0.357 

0.019 

0.611 

4 

0.007 

0.529 

0.440 

Table  6.  Example  scenario  probability  of  observation  from  optical  characteristic  observation 


Track 

Probability  of  RV 

Probability  of  Tank 

Probability  of  Decoy 

1 

0.784 

~0 

0.155 

2 

0.002 

0.703 

0.260 

3 

0.020 

0.350 

0.618 

4 

0.365 

0.018 

0.603 

Based  upon  these  distributions  and  the  measurements  collected  by  the  two  sensors,  EMs  were 
calculated  and  combined  with  Dempster’s  rule  of  combination.  The  resulting  PPs  are  contained 
in  Table  7.  Target  1  was  identified  as  having  the  greatest  PP  for  RV  and  was  selected  as  the 
target.  In  this  simulation,  target  1  corresponded  to  the  RV.  In  selecting  target  1,  the  system 
properly  discriminated  the  target  complex  and  selected  the  proper  object  for  engagement. 

Table  7.  Resulting  pignistic  probabilities  for  engagement 


Target 

PP(RV) 

PP(Tank) 

PP(Decoy) 

1 

0.868 

0.025 

0.107 

2 

0.028 

0.896 

0.076 

3 

0.123 

0.122 

0.755 

4 

0.149 

0.606 

0.193 

The  engagement  as  a  whole  was  successful  because  the  system  successfully  tracked  with  both 
sensors,  correlated  the  tracks,  and  properly  discriminated  the  complex  with  the  available 
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information.  While  this  engagement  was  successful,  the  system  could  still  fail  when  observing 
the  same  objects  and  produce  an  incorrect  discrimination  decision.  This  could  stem  from  the 
sensor  noise  corrupting  the  TOMs  or  random  variation  in  the  observed  characteristics. 
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5  Sensitivity  of  Sensor  Resolution  and  Characteristic  Overlap  on  Discrimination 


This  project  further  assessed  the  effect  of  varying  radar  angular  and  range  resolution, 
discrimination  characteristic  overlap  (the  amount  of  shared  area  between  the  probability  of 
observation  distributions  for  the  different  classes),  target  complex  size,  and  number  of  target 
complex  objects  on  the  performance  of  the  defense  system.  The  developed  simulation 
environment  was  utilized  to  perfonn  a  Monte  Carlo  simulation,  varying  the  parameters  listed 
above.  The  MATLAB  code  utilized  to  run  the  simulation  is  contained  in  Appendix  I  and  II.  The 
ultimate  measure  of  success  was  Pj.  Other  metrics  were  utilized  to  reflect  the  performance  of  the 
tracking,  correlation,  and  discrimination  elements  of  the  system  as  necessary  to  provide  a  more 
nuanced  perspective  on  the  performance  of  the  system. 

5.1  Effect  of  Radar  Resolution 

The  accuracy  of  the  radar  measurements  in  azimuth,  elevation,  and  range  can  affect  the 
radar’s  ability  to  produce  quality  tracks  for  the  engagement  of  the  threat  by  the  interceptor.  In 
order  to  assess  this  effect,  the  radar  range  and  angular  resolution  were  varied  by  +/-  10  percent  of 
the  baseline  values.  Elevation  resolution  was  held  at  2.3  times  azimuth  resolution  as  angular 
resolution  varied.  Range  and  angular  bias  were  set  with  a  uniformly  distributed  variation  with  a 
maximum  magnitude  equal  to  the  range  and  angular  resolution,  respectively.  Twenty- five 
different  combinations  of  values  were  tested  with  three  hundred  individual  trials  conducted  for 
each  individual  combination  of  values.  The  baseline  parameter  values  for  the  X-band  radar  are 
contained  in  Table  8. 


Table  8.  Baseline  radar  parameters 


Parameter 

Value 

Pulse  repetition  frequency  (Hz) 

45 

Azimuth  resolution  (deg) 

0.44 

Elevation  resolution  (deg) 

1.01 

Range  resolution  (m) 

10 

Figure  15  depicts  the  sensitivity  of  probability  of  discrimination  with  varying  range  and 
angular  resolution.  The  observed  probabilities  ranged  from  0.9  to  0.95.  However,  the  opposite  of 
the  anticipated  trend  was  noted,  as  Pd  actually  increased  as  range  and  angular  resolution 
degraded.  The  difference  between  minimum  and  maximum  Pd  varied  by  approximately  five 
percent,  a  relatively  small  variation  that  may  be  random  rather  than  reflective  of  system 
perfonnance.  In  addition,  the  implementation  of  the  Kalman  filter  and  the  correlation  process 
may  mitigate  the  effects  of  degrading  sensor  performance.  The  Kalman  filter  computed  a 
weighted  average  between  the  projected  position  of  a  track  and  the  observed  position.  By 
utilizing  a  projected  state,  the  filter  was  able  to  reduce  the  uncertainty  in  the  measurements  to 
less  than  the  uncertainty  associated  with  a  direct  measurement.  The  correlation  process  proved 
able  to  successfully  correlate  even  imperfect  TOMs.  Finally,  incorrect  correlations  between  RVs 
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and  decoys  may  have  reduced  the  impact  of  associating  the  wrong  tracks.  Since  the 
discrimination  characteristics  for  the  RV  and  decoys  are  relatively  similar,  a  correlation  between 
an  RV  and  a  decoy  may  still  produce  a  correct  discrimination  decision.  Since  the  complex 
contained  two  decoys,  one  RV,  and  one  tank,  it  is  more  likely  in  the  event  of  a  failed  correlation 
that  an  RV  and  decoy  track  be  paired  than  an  RV  and  tank. 
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Figure  15.  Variation  of  probability  of  discrimination  with  radar  range  and  angular  resolution 

Figure  16  depicts  the  variation  of  probability  of  correlation  with  radar  range  and  angular 
resolution.  The  observed  Pc  ranged  from  0.94  to  0.985.  With  small  exceptions,  Pc  increased  as 
angular  and  range  resolution  decreased.  This  follows  expectations,  as  improving  resolution 
should  improve  the  ability  of  the  radar  to  produce  higher  quality  tracks  that  will  more  closely 
match  the  scene  observed  by  the  KV  upon  target  acquisition. 
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Figure  16.  Variation  of  probability  of  correlation  with  radar  range  and  angular  resolution 

Figure  17  depicts  the  variation  of  mean  track  uncertainty  at  propagation  with  radar 
resolution.  Greater  track  uncertainty  at  propagation  can  lead  to  states  propagated  forward  to  the 
intercept  time  with  greater  uncertainty  and  errors,  leading  to  a  more  difficult  correlation  process. 
Degrading  angular  resolution  resulted  in  greater  track  uncertainty.  This  is  expected,  as  less 
precise  sensors  produce  less  certain  measurements.  Angular  resolution  had  a  much  larger  effect 
on  track  uncertainty  than  did  range  resolution.  Due  to  the  range  of  observation  (-2500  km),  a 
small  change  in  angular  resolution  produced  a  much  larger  effect  on  the  position  uncertainty.  For 
example,  to  produce  the  same  magnitude  in  Cartesian  position  uncertainty  as  a  0. 1  degree 
degradation  in  azimuth  resolution  at  a  range  of  2,500  km,  range  resolution  would  need  to 
degrade  by  4,360  m. 
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Figure  17.  Variation  of  mean  track  uncertainty  at  propagation  with  radar  resolution 

Figure  19  depicts  the  variation  in  Zero  Effort  Miss  distance  (ZEM)  with  radar  range  and 
angular  resolution.  ZEM  is  defined  as  the  Pythagorean  distance  between  the  aim  point  of  the  KV 
and  the  actual  location  of  the  selected  target  at  intercept  and  is  depicted  in  Figure  18. 


Figure  18.  Zero  Effort  Miss  distance  and  KV  field-of-view 
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Greater  ZEM  requires  a  more  maneuverable  KV  to  reach  the  target  and  a  greater  on-board  sensor 
FOV  to  detect  the  target.  As  angular  resolution  increased,  the  ZEM  increased.  Again,  as  shown 
in  Figure  19,  the  effect  of  angular  resolution  was  shown  to  be  much  greater  than  range 
resolution.  This  trend  in  ZEM  indicates  that  a  more  accurate  radar  will  place  smaller  demands  on 
the  KV  as  it  acquires  and  homes  on  the  target. 


Zero  Effort  Miss  Distance 


Angular  Resolution  (degrees) 


Figure  19.  Variation  of  zero  effort  miss  distance  with  radar  range  and  angular  resolution 

In  summary,  improving  radar  range  and  angular  resolution  improved  the  ability  of  the 
system  to  reduce  the  uncertainty  in  the  radar  tracks,  properly  correlate  the  radar  and  KV  TOMs, 
and  reduce  the  demands  on  KV  at  intercept.  However,  improved  radar  resolution  did  not 
translate  into  a  noted  improved  ability  to  discriminate  the  RV.  This  may  indicate  the  robust 
ability  of  the  Kalman  filter  and  the  correlation  algorithm  to  mitigate  the  effects  of  degrading 
sensor  performance.  Removing  the  Kalman  filter  from  the  simulation  and  propagating  the 
observed  states  from  the  radar  would  allow  for  direct  observation  of  the  effect  of  range  and 
angular  resolution  on  system  performance. 

5.2  Discrimination  Characteristic  Overlap  Area 

The  amount  of  overlap  between  the  probability  of  observation  distributions  for  different 
classes  is  indicative  of  the  likelihood  that  a  given  observation  could  indicate  different  classes.  An 
example  of  characteristic  overlap  is  shown  in  Figure  20. 
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Separation 


Class  Probability  of  Observation 


Figure  20.  Example  of  discrimination  characteristic  overlap  area  and  separation 

The  distributions  tend  to  move  together  for  two  reasons.  First,  less  capable  sensors  result  in 
distributions  with  higher  standard  deviations.  Greater  standard  deviations  tend  to  result  in  wider, 
more  overlapped  distributions.  Second,  objects  with  similar  characteristics  will  result  in 
distributions  with  more  closely  spaced  expected  values  and  tend  towards  more  overlap.  To  assess 
the  effect  of  the  amount  of  characteristic  overlap  area  on  Pd,  the  difference  between  the  expected 
values  of  the  classes  and  the  standard  deviation  of  the  measurements  were  varied  for  the  three 
discrimination  characteristics.  The  difference  between  the  means  was  varied  from  0. 1  to  5  and  the 
standard  deviations  were  varied  from  0.4  to  1.  The  amount  of  characteristic  overlap  was  calculated  for 
each  combination  of  mean  values  and  measurement  standard  deviation  with  1 000  trial  runs  for  each 
combination. 

Figure  2 1  depicts  the  variation  of  Pd  with  characteristic  overlap  area.  An  inverse  quadratic 
relationship  between  characteristic  overlap  area  and  Pd  was  observed,  with  an  r2  of  0.9836  for  the 
regression.  As  expected,  as  the  amount  of  overlap  area  increased,  the  ability  of  the  system  to  discriminate 
decreased.  Further,  this  relationship  was  strongly  quadratic,  indicating  that  Pd  decreases  more  rapidly  as 
the  overlap  area  increases.  Conversely,  decreasing  the  amount  of  overlap  area  will  yield  larger  increases 
in  Pd- 


Probability  of  Discrimination 
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Probability  of  Discrimination  v.  Characteristic  Overlap  Area 


Figure  21.  Variation  of  probability  of  discrimination  with  characteristic  overlap  area 
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6  Conclusions  and  Future  Research 


6.1  Conclusions 

Using  multiple,  dissimilar  sensors  operated  at  different  locations,  this  project  developed 
an  anchor  Target  Object  Map  (TOM)  with  least-squares  based  configuration  selection  for 
correlation  and  a  Dempster-Shafer  based  discrimination  algorithm  for  selecting  the  object  of 
interest  among  closely-spaced  objects  travelling  in  a  ballistic  trajectory.  A  simulation 
environment  was  developed  to  conduct  relevant  operational  scenarios  and  to  assess  the 
effectiveness  of  these  algorithms  with  varying  system  parameters.  This  project  found  that 
varying  +/-10%  of  the  baseline  radar  range  and  angular  resolution  capability  did  not  have  a 
significant  impact  on  the  probability  of  discrimination  (Pd)-  However,  a  0.1  degree  increase  in 
angular  resolution  resulted  in  a  6  percent  increase  in  the  zero  effort  miss  distance  that  increased 
the  required  KV  maneuverability  and  optical  sensor  suite  field  of  view  of  the  KV  at  intercept. 
This  project  also  found  a  strong  inverse  quadratic  relationship  between  Pci  and  discrimination 
characteristic  overlap,  which  is  defined  by  the  sensor  measurement  characteristics  for  different 
classes  of  objects.  In  order  to  achieve  a  Pd  of  0.9,  the  characteristic  overlap  area  could  not  exceed 
2.  Approaches  to  reduce  characteristic  overlap  area  should  be  explored  to  improve  the 
discrimination  performance  of  the  system. 

6.2  Future  Research 

The  simulation  capability  developed  in  this  project  could  be  improved  further  to  enhance 
the  performance  of  the  system  via  the  inclusion  of  an  Extended  Kalman  Filter22  for  tracking  and 
a  Mahalanobis  distance  correlation  algorithm.  In  addition,  the  fidelity  of  the  simulation 
environment  could  be  improved  by  modeling  the  KV  optical  sensor  utilizing  a  pinhole  optics 
model.  These  topics  have  been  explored  as  stretch  goals  of  the  project  but  were  not  able  to  be 
completed  in  a  timely  manner.  However,  the  ground  work  has  been  laid  and  can  be  employed  to 
further  enhance  the  multiple  sensor  discrimination  capabilities. 

6.2.1  Extended  Kalman  Filter  State  Estimation 

Due  to  the  nonlinear  nature  of  orbital  mechanics,  the  use  of  an  Extended  Kalman  Filter 
could  be  beneficial  in  improving  the  state  estimation  performance.  An  Extended  Kalman  Filter 
takes  the  previously  calculated  states  of  the  object  and  the  next  measurement  and  calculates  the 
true  states  of  the  object  at  the  given  time  step.  The  states  can  be  calculated  using  Equation 
(H).24 


xk  =  xk  +  Kk[zk-h(xk)]  (11) 

In  the  preceding  equation,  2fcare  the  estimated  states,  xk  are  the  previous  estimated  states, 
Kk  are  the  Kalman  gains  for  the  particular  measurement,  zi<  are  the  measured  states,  and  h(xk)) 
are  the  expected  measured  states  from  the  predicted  states.  xk  is  found  by  integrating  the 
nonlinear  equations  from  the  previous  states  to  the  next  state.  For  example,  Euler  integration  can 
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be  utilized  to  propagate  the  states  forward.  The  Kalman  gains  are  found  using  the  Riccatti 
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equations,  Equations  (  12  ),  (  13  ),  and  (  14  )  , 


Mk  =  <t>kPk^<t>Tk  +  Qk  (  12  ) 

Kk  =  MkHT(HMkHT  +  Rfc)-1  (  13  ) 

Pk  =  (/  -  Kk H)Mk  (  14  ) 


where  Mk  is  the  covariance  matrix  representing  the  errors  in  the  state  estimates  prior  to  an 
update,  and  Pk  is  the  covariance  matrix  following  an  update.  H  is  the  measurement  matrix 
relating  the  measurements  to  states.  Qk  is  the  process  noise  representing  uncertainty  in  the 
behavior  of  the  states,  and  Rk  is  the  discrete  measurement  noise  matrix  representing  the  variances 
of  each  measurement  source.  <f>k  is  the  fundamental  matrix,  found  via  the  Taylor-series  expansion 
given  in  Equation  (  15  ),  F  is  the  systems  dynamics  matrix  given  by  Equation  (  16  )  and  Ts  is  the 
sampling  time  between  measurements.  It  can  be  shown  that  the  higher  order  terms  can  be 
neglected  with  minimal  impact  on  the  accuracy  of  the  filter.  Since  the  system  is  nonlinear,  F  is  a 
first-order  approximation  to  linearize  the  equations. 
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6.2.2  Mahalanobis  Distance 

Due  to  sensor  uncertainty,  the  observed  location  of  each  object  may  differ  between  the 
two  sensors.  Mahalanobis  distance  can  be  used  to  calculate  the  correlation  between  the 
components  of  the  target  complex  and  is  computed  by  Equation  (  17  )  “  ,  where  r  is  the 
Mahalanobis  distance  from  the  feature  vector  x  to  the  mean  vector  im,  and  Cvis  the  covariance 
matrix  of  x. 


r2  —  {x_  —  m£)'Cx(x  —  m£)  (17) 

Mahalanobis  distance  is  utilized  as  it  removes  several  limitations  associated  with  using 
the  Euclidean  distance  mainly  that  it  accounts  for  the  scaling  of  the  coordinate  axes  and  corrects 
for  correlation  between  different  object  features.  The  Mahalanobis  distance  between  two  tracks 
is  denoted  by  rm-„,,  where  m  is  the  radar  track  number  and  n  is  the  optical  detection  track  number. 
The  selected  anchor  objects  will  always  be  correlated.  The  remaining  targets  will  be  correlated 
based  on  the  arrangement  that  minimizes  Er  for  the  TOM. 
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6.2.3  Pinhole  Optics  Model 

The  optical  sensor  onboard  the  KV  can  be  modeled  using  a  simple  pinhole  optics  model. 
The  pinhole  optics  model  involves  defining  a  detector  with  a  certain  number  of  pixels  and 
projecting  this  sensor  plane  through  a  “pinhole”  to  measure  the  range  and  bearing  of  the  objects 
in  the  field  of  view  (FOV).  Each  pixel  can  be  assumed  to  have  square  dimensions.  The  sensor 
FOV  was  then  divided  by  the  number  of  pixels,  known  as  the  instantaneous  FOV  (iFOV)  for 
each  pixel.  The  FOV  and  iFOV  of  the  sensor  array  will  be  varied  to  assess  their  impact  on 
system  performance. 
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Appendix  I:  Monte  Carlo  Simulation  Coding 


MIDN  Samuel  S.  Lacinski,  25  April  2014 

Project  Abstract 

This  project  examines  the  selection  of  an  object  of  interest  (01)  from  a  cluster  of  closely-space 
objects  (CSOs)  utilizing  the  returns  from  multiple  sensors. 


%  This  project  is  conducted  in  the  SI  system  of  units.  Care  is  taken 
%  between  sensor  measurements,  taken  in  units  of  meters,  and  orbital 
%  mechanics  calculations,  in  units  of  kilometers. 

%  This  program  has  five  principle  components:  1.  Turth  Target  Complex 
%  Generation,  2.  Radar  Observation,  3.  Optical  Observation, 

%  4.  Target-Object  Map  (TOM)  Correlation,  and  5.  Discrimiation  Algorithm. 
%  In  element  will  be  explained  in  detail  in  the  respective  section  of 
%  coding. 


Workspace  Preparation 


clear 

cl  c 

format  compact 

%  Declare  Path 

global  currPath; 
currPath  =  pwd ; 

path([currPath 

'/CoordinateTransforms/'] ,path) 

path([currPath 

'/Discrimination/']  ,path) 

path([currPath 

'/Filter/']  .path) 

path([currPath 

'/Sensors/'] ,path) 

path([currPath 

'/TimeAndAngles/'] .path) 

path([currPath 

'/Trajectory/'] ,path) 

path([currPath 

'/DataFiles/'] .path) 

path([currPath 

'/EClFilter/']  .path) 

path([currPath 

'/SEZFilter/']  .path) 

path([currPath 

' /AddAxi s/'] , path) 

global  plotk; 

plotk  =  []  ; 

%  plot  i ndex 

Compute  Truth  Trajectories  and  Access  Data 

%  Earth  Paratemeters 
global  mu  R_e  w; 
mu  =  3 . 986004418e5 ; 

R_e  =  6371; 


%  kmA3/sA2,  Earth  Gravitational  Parameter 
%  km,  Earth  radius 
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w  =  [0;0;7.292115856e-5]  ;  %  rad/s,  Earth  angular  rotation 


%  Target  Complex  Components 
num_l  =1;  %  n 


%  number  of  re-entry  vehicles 
%  number  of  tanks 
%  number  of  spherical  decoys 
%  number  of  conical  decoys 
%  number  of  debris  objects 


num_2  =  1 
num_3  =  2 
num_4  =  0 
num_5  =  0 


target_comp  =  [num_l  num_2  num_3  num_4  num_5] ; 

%  Target  components  vector 
num_tar  =  sum(target_comp) ;  %  number  of  targets 

%  Define  target  parameters 
global  num_class; 

num_class  =3;  %  number  of  object  classes 

v_const  =  .01; 

%  global  UTC_launch  JD_launch  file; 

%  file  =  '  PositionTrackData2.xlsx'  ;  %  data  file 

%  Define  conops 
global  int_aq; 

int_aq  =  500;  %  km,  range  of  acquisition 

%  Define  Xl  position 

radar_lat  =  29.6837;  %  deg,  radar  latitude 

radar_long  =  -165.192;  %  deg,  radar  longitude 

radar_alt  =  0;  %  km,  radar  altitude 

Xl.data  =  [radar_lat  radar_long  radar_alt]; 

Truth  Target  Complex 

This  section  of  code  generates  the  ballistic  trajectory  of  the  components  of  the  CSOs.  This 
project  is  focused  on  the  free-flight  portion  of  the  CSOs.  This  section  of  code  also  generates  the 
true  object  charateristics.  Each  object  has  six  total  characteristics:  a,  b,  c,  alpha,  beta,  and 
gamma,  a,  b,  and  c  are  radar  observable  charateristics  while  alpha,  beta,  and  gamma  are  optical 
observable  charateristics. 

The  generated  information  is  contained  in  a  structured  array,  target  contains  a  Nxl  cell  matrix, 
where  N  is  the  true  number  of  target  elements.  The  trajectory  element  of  each  cell  contains  the 
julian  day  and  corresponding  J2000  position  and  velocity  elements  of  the  target  elements. 

%  Compute  trajectory  of  target  elements 

[complex, connect]  =  Scenariolnitial (target_comp,Xl.data,v_const) ; 
app  =  connect. app; 
root  =  connect . root ; 
scenario  =  connect. scenario; 


42 


%  Calculate  discrimination  characteristics  and  obtain  epoch 
in  =  1;  %  target  index 

for  k  =  l:num_class 

for  n  =  1: target_comp(k) 
if  k  ==  1 

complex{in}.char  =  [14]; 
el  seif  k  ==  2 

complex{in}.char  =[3  6]; 

else 

complex{in}.char  =  [2  5]; 

end 

in  =  i n+1 ; 

end 

end 

%  Declare  expected  observation  for  target  classes 
global  mean_char; 
mean_char  =  [1  3  2;4  6  5]; 

%  Determine  initial  launch  time 
global  UTC_launch  JD_launch; 

UTC_launch  =  complex{l} . UTClaunch ; 

JD_launch  =  juliandate(UTC_launch) ; 

%  Import  Access  Data 

Xl . access  =  cel  1  (num_tar ,  1)  ;  %  initialize  access  cell  matrix 

global  JD_hitch 

in  =  1;  %  reset  index 

radar  =  connect. Xlstk;  %  select  radar  facility 
for  k  =  l:num_tar  %  Loop  over  all  complex  objects 
%  Compute  access  data  for  each  object 
if  k  ==  1 

Tar  =  connect. RVs; 
el  seif  k  ==  2 

Tar  =  connect. Tanks; 
el  seif  k  ==  3 

Tar  =  connect. Spheres; 
el  seif  k  ==  4 

Tar  =  connect. Cones; 

el  se 

Tar  =  connect. Debrs; 

end 

for  n  =  l:size(Tar,l) 

satid  =  Tar{n,l};  %  object  identifier 

sat  =  Tar{n,2};  %  object 

%  Compute  Access 

access  =  radar. GetAccessToOb ject(sat) ; 
access . ComputeAccess ; 

accessDP  =  access. DataProviders.ltemC'AER 
Data ' ) . G  roup . Item( ' Def aul t ' ) . Exec (scenari o . startti me , scenari o . stopti me , 1) ; 
T  =  accessDP. DataSets. GetDataSetByName('Time') .Getvalues; 
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T  =  UTC2Scenario(T)  ; 

T  =  jul i andate(T) ; 

az  =  cel  1 2mat(accessDP. DataSets .GetDataSetByName( 'Azimuth ' ) .Getval  ues)  ; 

el  =  cel  12mat(accessDP.  DataSets.  GetDataSetByNameC  Elevation')  .Getval  ues)  ; 

r  =  cell2mat(accessDP. DataSets. GetDataSetByNameC ' Range ' ) .Getvalues) ; 

azdot  =  cell2mat(accessDP. DataSets. GetDataSetByNameC 'AzimuthRate ' ) .Getvalues) ; 

eldot  =  cell2mat(accessDP. DataSets. GetDataSetByName('ElevationRate') .Getvalues) ; 

rdot  =  cell2mat(accessDP. DataSets. GetDataSetByNameC ' RangeRate ' ) .Getvalues) ; 

%  Save  data 

Xl. access{i n}  =  [T  az  el  r  azdot  eldot  rdot]; 

%  Find  di sconti nuni ty 

DD_hitch(in,  :)  =  FindAERDiscon(Xl.access{in}(:  ,1)  ,xl.access{in}(:  ,2))  ; 
in  =  in+1;  %  update  index 

end 

end 

%  Close  STK 
app. Quit()  ; 

%  Load  scenario  values 
load ('Scenario. mat') 


Declare  Range  of  Variables  for  Iteration 

Default  Values:  sigma_res  =  .0458  deg  sigma_rng  =  10/1  e3  km 


res  =  0.44;  %  deg,  baseline  angular  resolution 

rng  =  10/le3;  %  km,  baseline  range  resolution 

sigma_res  =  [.9*res  .95*res  res  1.05*res  1.1* res];  %  deg,  angular  resolution  to  be  tested 

sigma_rng  =  1.05*rng;  %  km,  range  resolution  to  be  tested 


Perform  Iteration 

This  section  of  code  implements  a  Monte  Carlo  simulation  to  assess  the  effect  of  the  variation  of 
the  parameters  of  interest. 


%  Save  globals  into  structured  variable 

globals.plotk  =  plotk; 

globals. mu  =  mu; 

global s.R_e  =  R_e; 

globals. w  =  w; 

globals.num_class  =  num_class; 
global  s.int_aq  =  int_aq; 
global s.mean_char  =  mean_char; 
global s.UTC_launch  =  UTC_launch; 
global s.DD_launch  =  JD_launch; 
global s . DD_hitch  =  JD_hitch; 


loops  =  300; 


%  number  of  iterations  at  each  combination 
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%load( ' LastRun .mat ' ) 
n_last  =  1; 
m_last  =  1; 

results  =  zeros(loops,10)  ; 

P  =  zeros(length(sigma_res) ,10,length(sigma_rng)) ; 

for  n  =  n_last: length (si gma_rng) 
if  n  ==  n_last 

m_l oop  =  m_l ast ; 

el  se 

m_loop  =  1; 

end 

for  m  =  m_loop:length(sigma_res) 
m 
n 

parfor  k  =  1: loops 
test  =  0; 

iter  =0;  %  number  of  attempts 

while  test  ==  0  &&  iter  <  10 
try 

%  Run  simulation 

results(k,  :)  =  CSOMultiSensorDiscrimLoop(complex,Xl,globals,  .  .  . 

sigma_res(m) , si gma_rng(n)) ; 
test  =  1; 
catch 

%  loop  repeats 

dispC'Error  detected.  Repeating  data  point.') 
iter  =  iter+1; 

end 

end 

if  iter  >=  10 

results(k,:)  =  NaN(l,8); 

end 

end 

P(m,:,n)  =  sum(results)/loops ; 
save ('LastRun. mat' ,'m','n','P') 
end 

end 


Published  with  MATLAB®  R2014a 
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Appendix  II:  Simulation  Master  Code 


function  [results]  =  CSOMultiSensorDiscrimLoop(complex,Xl,globals,sigma_res,sigma_rng) 


Declare  Variables 

global  plotk  mu  R_e  w  num_class  i nt_aq  mean_char  UTC_launch  DD_launch  JD_hitch 
plotk  =  global  s.  plotk; 
mu  =  global  s.  mu; 

R_e  =  global s.R_e; 
w  =  globals.w; 

num_class  =  globals.num_class; 
int_aq  =  global s.int_aq; 
mean_char  =  global s.mean_char; 

UTC_launch  =  global s.UTC_launch; 

JD_launch  =  gl obal s . JD_1 aunch ; 

JD_hitch  =  global s.JD_hitch; 

Error  using  CSOMultiSensorDiscrimLoop  (line  4) 

Not  enough  input  arguments. 

Create  STK  Scenario 

%  Create  scenario 

app  =  actxserverC'STKlO. Application') ;  %  open  STK 
root  =  app. personality^; 

scenario  =  root. Children. New('eScenario' , 'Kwaj2seattle') ; 

%  Set  start  and  stop  time 

scenario. SetTimePeriod('l  Sep  2014  00:00:00','l  Sep  2014  01:00:00'); 

%  Save  to  connect  variable 
connect. app  =  app; 
connect. root  =  root; 
connect . scenario  =  scenario; 


Radar  Observation 

This  section  of  code  simulates  the  observation  of  the  target  complex  by  a  radar  system.  The  radar 
observes  the  target  for  a  given  period  and  then  generates  a  TOM  that  is  projected  to  the  "eyes- 
open"  time  for  the  optical  sensor.  The  return  from  the  complex  is  the  target  position,  velocity, 
and  the  radar  signal  return  from  the  target. 

This  project  assumed  that  the  radar  would  add  a  random  error  to  error  element's  truth  location. 
This  error  is  based  on  the  radar  system  parameters.  In  addition,  the  radar  signal  return  will  have 
some  random  noise  added  to  it  based  upon  the  internal  noise  of  the  system. 
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The  radar  simulated  in  this  scenario  is  a  nominal  X-Band  midcourse  tracking  and  discrimination 
radar,  termed  the  XI.  The  radar  tracking  parameters  are  stored  in  a  structured  array  that  contains 
the  parameters  of  the  system  in  .parameters  and  target  tracks  in  .tracks. 


%  Define  the  radar  system 

parameters 

Xl.sigma_rng  =  sigma_rng; 

%  km,  range  standard  deviation 

Xl.ang_res  =  sigma_res; 

%  deg,  angular  resolution 

Xl.sigma_a  =0.5; 

%  characteristic  a  standard  deviation 

Xl. PRF  =  45 ; 

%  Hz,  pulse  repetition  frequency 

XI. Pd  =  1; 

%  Probability  of  detection 

Xl.lat  =  Xl.data(l); 

%  deg,  radar  latitude 

Xl.long  =  Xl.data(2); 

%  deg,  radar  longitude 

Xl.alt  =  Xl.data(B); 

%  km,  radar  altitude 

%  Generate  target  tracks 

[Xl,KV, complex, connect]  = 
dispC'Radar  track') 

radarOBVwSTK(Xl, compl  ex , connect)  ; 

Optical  Observation 

This  section  of  code  simulates  the  observation  of  the  target  complex  by  a  space-based  optical 
system.  The  sensor  platform  is  also  moving  on  a  ballistic  trajectory  and  observes  the  complex  for 
"eyes  open"  time  to  the  decision  point.  At  the  decision  point,  the  observations  are  used  to 
generate  a  TOM.  The  sensor  observes  the  2D  position  of  the  CSOs  in  the  focal  plane  as  well  as 
characteristics  of  interest  for  discrimination  purposes. 

This  project  assumes  the  optical  sensor  would  add  a  random  error  to  each  element's  truth 
location,  based  on  the  system  parameters.  In  addition,  the  observed  characteristics  would6  have 
an  error  resulting  from  the  internal  noise  of  the  system. 

%  Define  system  parameters 
KV. si gma_l  =0.5; 

%  Observe  target  and  generate  TOM 
kv  =  kvobv_stk(x1,kv, complex) ; 
dispC'KV  Acquired') 


TOM  Correlation 

This  section  of  code  correlates  the  final  TOMs  generated  by  the  radar  and  optical  sensors.  The 
correlation  parameter  generated  is  based  upon  the  "2-D  correlation"  of  the  coordinates  and 
measured  velocities  of  the  objects. 

This  project  assumes  no  mismatches. 
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[optical_obv,TOM_data]  =  TOMCorr(Xl,KV,complex) ; 
if  ~i semptyCpl otk) 

disp(['OV  selected  TOM  configuration:  1 ,num2str(TOM_data(2))]) 
disp(['r  =  ' ,num2str(TOM_data(l) ,4)]) 


Discrimination  Algorithm 

This  section  of  code  fuses  the  discrimination  information  from  the  two  sensors  and  produces  a 
discrimination  decision.  The  algorithm  fuses  the  sensor  inputs  based  upon  the  TOM  correlation 
process  in  the  previous  section  of  code.  Objects  selected  as  matches  by  the  correlation  process 
will  have  the  sensor  returns  matched  in  the  object  data  base. 


%  Combine  observations  based  on  correlation  process. 
char_data  =  [Xl.char  KV. char( : , 2 : end)] ; 

%  Assign  evidence  measure  based  on  each  observation 

num_obj  =  size(char_data,l) ;  %  number  of  correlated  objects 

num_char  =2;  %  number  of  characteristics  measured 

char_sigma  =  [Xl.sigma_a  KV.sigma_l]; 

PPs  =  DempsterShaferOpen(mean_char,char_sigma,char_data) ; 

%  Identify  Object  of  Interest 

[F_imax,h]  =  max(PPs(:  ,1))  ; 
if  ~i  semptyCpl  otk) 

di sp( [ 'The  OI  is  ' , num2str(h)]) 

di sp( [' Interest  factor  is  1 , num2str(F_imax)] ) 

end 


Save  Results 

This  section  of  coding  determines  if  correlation,  discrimination,  and  intercept  were  successful 
and  saves  information  of  interest. 


%  Determine  if  discrimination  was  successful  (value  of  1  indicates  success) 
if  h  ==  1  &&  ~isnan(TOM_data(l)) 

Disc_suc  =  1; 

el  se 

Disc_suc  =  0; 

end 

%  Determine  if  correct  anchors  were  selected 
if  Xl.anchor_id  ==  2  &&  ~isnan(TOM_data(l)) 
radar_anc  =  1; 

el  se 

radar_anc  =  0; 


end 

if  KV.anchor_id  ==  2  &&  ~i snan(TOM_data(l)) 
kv_anc  =  1; 

el  se 

kv_anc  =  0 ; 

end 

%  Determine  if  KV  selected  correct  target  independently 
if  KV. target  ==  1  &&  ~i snan(TOM_data(l)) 

KV_int  =  1; 

el  se 

KV_int  =  0; 

end 

%  Determine  if  intercept  was  successful 
if  strcmp(complex{h,l}.lD, 'RV')  &&  ~isnan(TOM_data(l)) 
lnt_suc  =  1; 

el  se 

lnt_suc  =  0; 

TOM_data(l)  =  -1; 

end 

%  Save  results  (format  -  [Corr  Di scrim  Int  Corr. coefficient  PP_RV]) 
results  =  [TOM_data(3)  Disc_suc  lnt_suc  TOM_data(l)  F_imax  radar_anc 
mean(Xl.TU_min)  KV_int  KV.zem]; 
di sp( ' Intercepted ' ) 
disp('  ') 

%  Close  STK 
if  isempty(plotk) 
app. Quit()  ; 

end 


Published  with  MATLAB®  R2014a 
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Appendix  III:  Kalman  Filter  Code 


function  [xkl, Pkl,Kk,TU,TU_v]  =  KalmanSEZ(xO, P0,Rk,zk,dt) 

%KalmanSEZ  performs  one  step  of  a  linear  Kalman  filter  on  states  in  the  SEZ 
%f rame . 

%  The  states  are  expected  in  the  form  [S,E,Z,Sdot,Edot,Zdot] ' ,  where 
%  position  is  in  km  and  velocity  is  in  km/s.  Pk  is  a  6x6  covariance 

%  matrix  in  kmA2.  Rk  is  the  measurement  covariance  matrix  also  in  kmA2. 

%  Zk  is  the  measurement  matrix,  expected  as  [S,E,Z],  as  velocity  is 
%  assumed  to  not  be  measured.  If  Zk  is  empty,  the  filter  is  propagated 
%  without  a  measurement  vector.  T  is  the  propagation  time  in  seconds. 

global  mu  R_e; 

Rk  =  [Rk(l, 1)  0  0; . . . 

0  Rk(2 , 2)  0; 

0  0  Rk(3,3)]  ; 

%  State  transition  matrix 
phi  =  eye (6, 6)  ; 
phi (1,2)  =  dt; 
phi (3,4)  =  dt; 
phi (5,6)  =  dt; 

%  Input  matrix 
G  =  [dtA2/2  0  0; . . . 
dt  0  0; . . . 

0  dtA2/2  0;  .  .  . 

0  dt  0; . . . 

0  0  dtA2/2; . . . 

0  0  dt]  ; 

%  Process  matrix 

q  =  le-6;  %  Noise  spectral  density 

Q  =  q* [dtA3/3  dtA2/2  0  0  0  0;  .  .  . 
dtA2/2  dt  0  0  0  0; . .  . 

0  0  dtA3/3  dtA2/2  0  0; . . . 

0  0  dtA2/2  dt  0  0;  .  .  . 

0000  dtA3/3  dtA2/2; . . . 

0  0  0  0  dtA2/2  dt]  ; 

%  Measurement  Matrix 
H  =  zeros(3 , 6) ; 

h (1 , 1)  =  1; 

H (2 , 3)  =  1; 

H (3 , 5)  =  1; 

%  Gravitational  Acceleration  Matrix 

g  =  -mu/norm([x0(l)  X0(3)  X0(5)+R_e] )A2 ;  %  estimated  gravitational  constant 

u_k  =  [0;0;g] ; 
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%  Best  Estimate  States 

Xk  =  phi*xO+G*u_k;  %  km,  km/s,  best  estimate  states 

Pk  =  phi *pO*phi '+Q;  %  kmA2,  best  estimate  covariance 

%  Update  Estimates  -  omit  if  no  measurement  included 
if  ~isempty(zk) 

Kk  =  Pk*H'*inv(H*Pk*H'+Rk)  ;  %  Kalman  gain  matrix 
Xkl  =  Xk+Kk* (zk-H*Xk) ;  %  km,  km/s,  updated  estimates 

Pkl  =  (eye(6,6)-Kk*H)*Pk;  %  kmA2 ,  updated  covariance  matrix 

else  %  Use  propagated  estimates  as  updates 
Kk  =  NaN(6 , B)  ; 
xkl  =  xk; 

Pkl  =  Pk; 

end 

%  Calculate  Track  Uncertainty  at  measurement  instance 
TU  =  norm([Pkl(l,l)  Pkl(2,2)  Pkl(3,3)]); 

TU_v  =  norm([Pkl(4,4)  Pkl(5,5)  Pkl(6,6)]); 

end 


Published  with  MATLAB®  R2014a 
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Appendix  IV:  MATLAB-STK  Interface  for  Calculating  Interceptor  Trajectory 


function  [radar , KV, connect]  =  lntercept(radar,intercept_data, connect) 

%lntercept  This  function  determines  the  time  of  intercept  for  a  given 

%threat  trajectory. 

%  Output  is  in  JD. 


Variable  Definiation 

global  int_aq; 

num_tar  =  size(radar.trackSEZ,l) ;  %  number  of  active  tracks 

%  STK  scenario  data 
app  =  connect. app; 
root  =  connect . root ; 
scenario  =  connect. scenario; 

Error  using  Intercept  (line  10) 

Not  enough  input  arguments. 

Select  State  Estimates 

State  estimates  with  the  smallest  average  covariance  nonnal  for  all  tracked  objects  are  utlized  for 
determining  the  object  orbit. 


%  %  Determine  number  of  loops 
%  loops  =  inf; 

%  for  k  =  l:size(radar.trackSEZ,l) 

%  det  =  size(radar.trackSEZ{k,2},3) ; 

%  if  det  <  loops 

%  loops  =  det; 

%  end 

%  end 

%  ave  =  zeros(loops,l) ; 

% 

%  %  Loop  over  each  track  indicidence 
%  for  k  =  l:loops 
%%  disp(['Loop:  ' ,num2str(k)]) 

%  aves  =  zeros (num_tar, 1) ; 

%  %  Average  covariance  for  each  track 

%  for  n  =  l:num_tar 

%%  disp(['Target:  ' , num2str(n)]) 

%  aves(n,l)  =  norm(radar.trackSEZ{n,2}(: , : ,k)) ; 

%  end 

%  ave(k,l)  =  norm(aves);  %  kmA2 ,  average  covariance  for  all  tracks 
%  end 
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%  Select  point  of  minimum  track  uncertainty 

h  =  zeros(l, num_tar) ; 

radar. TU_min  =  zeros (num_tar, 1) ; 

for  k  =  l:num_tar 

[ radar. TU_mi n (k) ,h(k)]  =  mi n(radar . trackSEZ{k, 4}( : , 3)) ; 

end 

%  Select  states  for  orbit  determination 

%[~, h]  =  mi n(ave) ; 

state_prop  =  zeros(num_tar,8) ; 

for  k  =  l:num_tar 

statesSEZ  =  radar. trackSEZ{k,l}(h(k) , :) ;  %  states  for  propagation 
%  Convert  states  to  ECI 

3D  =  statesSEZ(l,2) ;  %  Julian  day  for  state 

theta_LST  =  LST(JD, radar. long) ;  %  deg,  local  sidereal  angle 

[r_ECl, v_ECl]  =  SEZ2 ECI (statesSEZ (3 : 5) ' , statesSEZ(6:8) ' , radar .1  at , . . . 

theta_LST);  %  km,  km/s ,  ECI  states 

state_prop(k, :)  =  [statesSEZ(l: 2)  r_ECl '  v_ECl']; 

end 


Project  Trajectories 

Using  the  selected  state  estimates,  the  tracks  are  projected  forward  until  impact  to  detennine  the 
trajectory  of  the  targets  after  tracking. 


tars  =  cell  (num_tar , 2)  ;  %  initialize  object  handle  matrix 

root . ExecuteCommand( ' Setuni ts  /  KM  SEC  UTCG'); 
for  k  =  l:num_tar 

tarid  =  ['TAR' ,num2str(k)] ;  %  target  label 

tar  =  scenario. Children. New(' esatellite' , tarid) ;  %  target  STK  object 

%  Initialize  target  with  state  estimates 

[year, month, day, hour, minu, sec, ~,~]  =  jul i an2greg(state_prop(k, 2)) ;  %  UTC  time 

%  Generate  time  string 
T  =  [year  month  day  hour  minu  sec]; 
radar. prop_start(k,l)  =  juliandate(T) ; 

T  =  UTC2str(T) ; 

%  Add  object 

cmd  =  ['setstate  */Satel 1 i te/ ' , tari d ,  '  Cartesian  TwoBody  ' ,T, '  scenario . StopTime , 
5  J2000  ' ,T, 1  '  ,num2str(state_prop(k,3)) , '  '  ,num2str(state_prop(k,4)) , . . . 

'  ' , num2str(state_prop(k, 5)) , '  ' ,num2str(state_prop(k,6)) , . . . 

'  ' , num2str(state_prop(k, 7)) , '  ' ,num2str(state_prop(k,8))] ; 
root . ExecuteCommand (cmd) ; 

tars{k,l}  =  tarid;  %  save  object  name 

tars{k,2}  =  tar;  %  save  object  handle 

%  Extract  projected  positions 
satDP  =  tar. DataProvi ders . Item( ' Cartesian 
Posi  ti  on ' ) .Group . Item( ' J2000 ' ) . Exec(scenari  o . startti  me , scenari  o . stopti  me , 1) ; 

T  =  satDP. DataSets. GetDataSetByName('Time') .Getvalues; 

T  =  UTC2scenario(T)  ; 

T  =  jul i andate(T) ; 
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x  =  cell2mat(satDP. DataSets. GetDataSetByName('x') .Getvalues) ; 
y  =  cell2mat(satDP. DataSets. GetDataSetByName('y') .Getvalues) ; 
z  =  cell2mat(satDP. DataSets. GetDataSetByName('z') .Getvalues)  ; 
satDP  =  tar. DataProviders.ltemC' Cartesian 
Velocity') . G roup. Item ( ' J 2000') .Exec(scenario.starttime,scenario.stoptime,l) ; 
xdot  =  cell 2mat(satDP . DataSets. GetDataSetByName('x') .Getvalues) ; 
ydot  =  cell2mat(satDP. DataSets. GetDataSetByName('y') .Getvalues)  ; 
zdot  =  cell 2mat(satDP . DataSets. GetDataSetByName('z') .Getvalues) ; 

%  Save  data 

radar. prop_traj{k,l}  =  [T, x,y,z, xdot, ydot, zdot] ; 

end 


Determine  Intercept  Point 

Based  on  the  provided  altitude  in  intercept  data,  the  function  will  determine  the  time  and 
location  of  intercept. 


int_alt  =  intercept_data. altitude;  %  Desired  intercept  altitude 

%  Identify  earliest  time  object  breaks  intercept  altitude 
T_int  =  inf;  %  initialize  intercept  time 

for  k  =  l:num_tar 
%  Select  object 
tarid  =  tars{k,l}; 
tar  =  tars{k,2} ; 

%  Compute  altitude  over  flight 
tarDP  =  tar.DataProviders.ltem('LLA 

State ' ) .Group. Item( ' Fixed ' ) . Execfscenari  o . starttime , scenario. stopti  me , 1) ; 
T  =  tarDP. DataSets. GetDataSetByName('Time') .Getvalues; 

T  =  UTC2scenario(T) ;  %  UTC 

T  =  jul i andate(T) ;  %  Julian  Day 

alt  =  cell2mat(tarDP.DataSet.GetDataSetByName('Alt') .Getvalues) ; 

%  Compute  time  at  intercept  altitude 
T_1 ow  =  interpl(alt,T,int_alt) ;  %  Julian  Day 
%  Update  intercept  time,  if  necessary 
if  T_1  ow  <  T_int 
T_int  =  T_low; 

end 


inter_rv  =  zeros(num_tar,7) ; 

%  Determine  cartesion  position  and  velocity  at  intercept  time 
for  k  =  l:num_tar 
%  Select  target 
tarid  =  tars{k,l}; 
tar  =  tars{k,2} ; 

%  Obtain  Cartesian  position  and  velocity 
satDP  =  tar. DataProviders.ltemC' Cartesian 
Position') .Group. Item( ' J2000 ' ) . Exec(scenario . starttime, scenario. stoptime, 1) ; 
T  =  satDP. DataSets. GetDataSetByName('Time') .Getvalues; 

T  =  UTC2scenario(T)  ; 

T  =  jul i andate(T) ; 

x  =  cell2mat(satDP. DataSets. GetDataSetByName('x') .Getvalues)  ; 
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y  =  cell2mat(satDP. DataSets. GetDataSetByName('y') .Getvalues) ; 
z  =  cell2mat(satDP. DataSets. GetDataSetByName('z') .Getvalues) ; 
satDP  =  tar. DataProviders.ltemC' Cartesian 
Velocity') . Group. Item ( ' J 2000') .Exec(scenario.starttime,scenario.stoptime,l)  ; 
xdot  =  cell 2mat(satDP . DataSets. GetDataSetByName('x') .Getvalues) ; 
ydot  =  cell 2mat(satDP . DataSets. GetDataSetByName('y') .Getvalues) ; 
zdot  =  cell2mat(satDP. DataSets. GetDataSetByName('z') .Getvalues) ; 

%  Interpolate  position  and  velocity  at  intercept  time 
X  =  interpl(T,x,T_int) ; 

Y  =  interpl(T,y,T_int) ; 

Z  =  interpl(T,z,T_int) ; 

Xdot  =  interpl(T,xdot,T_int) ; 

Ydot  =  interpl(T,ydot,T_int) ; 

Zdot  =  i nterpl(T, zdot ,T_i nt) ; 

%  Save  data 

inter_rv(k, :)  =  [T_int  X  Y  Z  Xdot  Ydot  Zdot]; 

end 

%  Determine  aim  point  (aim  point  taken  as  average  of  all  objects  position 
%  at  intercept  time) 

aim_point  =  [T_int  mean(i  nter_rv(  :  , 2))  mean(i nter_rv( :  ,  3))  mean(inter_rv(:  ,4))]  ; 
KV.aim_point  =  aim_point; 


Determine  Interceptor  Trajectory 

This  segment  of  code  calculates  the  trajectory  necessary  for  the  interceptor  to  converge  on  the 
target. 


%  Create  interceptor  object 

KVl  =  scenario. Chi Idren.NewCeSatellite'  ,  'KVl')  ; 

%  Create  target  sequence  and  remove  default  initial  state  and  propagat 
root . ExecuteCommand( ' Astrogator  */Satel 1 i te/KVl  Insertsegment 
Mai nSequence . SegmentLi st . Ini ti al_State  Target_Sequence 1 ) ; 
root . ExecuteCommand( 'Astrogator  */Satel 1 i te/KVl  DeleteSegment 
Mai nSequence . SegmentLi st . Ini ti al_State  Target_Sequence ' ) ; 

root . ExecuteCommand( 'Astrogator  */Satel 1 i te/KVl  DeleteSegment  Mai nSequence . SegmentLi st . Propagate 
Target_Sequence ' ) ; 

%  Create  launch,  manuever,  and  propagate  sequence 

root . ExecuteCommand( 'Astrogator  */Satel 1 i te/KVl  Insertsegment 

Mai nSequence . SegmentLi st .Target_Sequence. SegmentLi  st . -  Propagate ' ) ; 

root . ExecuteCommand( 'Astrogator  */Satel 1 i te/KVl  Insertsegment 

Mai nSequence . SegmentLi st .Target_Sequence. SegmentLi st . Propagate  Launch ' ) ; 

root . ExecuteCommand( 'Astrogator  */Satel 1 i te/KVl  Insertsegment 

Mai nSequence . SegmentLi st .Target_Sequence. SegmentLi st . Propagate  Maneuver ' ) ; 

%  Set  constants  for  launch 

root . ExecuteCommand( 'Astrogator  */Satel 1 i te/KVl  Setvalue 

Mai nSequence . SegmentLi st .Target_Sequence . SegmentLi st . Launch .Ti meof Fl i ght  63  sec ' ) ; 
root . ExecuteCommand( 'Astrogator  */Satell i te/KVl  Setvalue 

Mai nSequence . SegmentLi st .Target_Sequence. SegmentLi st . Launch . Burnout. Fi xedveloci ty  8 . 3  km/sec ' ) ; 
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root . ExecuteCommand( ' Astrogator  */Satell i te/KVl  Setvalue 

Mai nsequence . Segment Li st .Target_Sequence . Segment Li st . Launch . Launch . Geodeti c . Lati tude  B4 . 7561 
deg')  ; 

root . ExecuteCommand( 'Astrogator  */Satell i te/KVl  Setvalue 

Mai  nsequence .  Segment Li  st .Target_Sequence . Segment Li st . Launch . Launch . Geodeti c . Longi tude  -120.626 
deg')  ; 

root . ExecuteCommand( 'Astrogator  */Satell i te/KVl  Setvalue 

Mai  nsequence.  Segment Li  st.Target_Sequence. Segment List. Launch. Launch. Geodetic.  Altitude  0.01  km')  ; 
root . ExecuteCommand( 'Astrogator  */Satell i te/KVl  Setvalue 

Mai nsequence . Segment Li st .Target_Sequence . Segment Li st . Maneuver . Impul si veMnvr . Atti tudecontrol 
Thrust  Vector' ) ; 

%  Configure  target  sequence  control  parameters 

root . ExecuteCommand( 'Astrogator  */Satel 1 i te/KVl  AddMCSSegmentControl 

Mai nsequence. Segment Li st.Target_Sequence. Segment List. Launch  Burnout. LaunchAzDRDAlt. Altitude') ; 

root . ExecuteCommand( 'Astrogator  */Satel 1 i te/KVl  SetMCSControl Val ue 

Mai nsequence . SegmentLi st .Target_Sequence . Profi 1 es . Di fferenti al_Corrector  Launch 

Burnout. LaunchAzDRDAlt. Altitude  Active  true'); 

root . ExecuteCommand( 'Astrogator  */Satell i te/KVl  Setvalue 

Mai nsequence. Segment Li st.Target_Sequence. Segment List. Launch. Burnout. LaunchAzDRDAlt. Altitude  225 
km') ; 

root . ExecuteCommand( 'Astrogator  */Satel 1 i te/KVl  AddMCSSegmentControl 

Mai nsequence . SegmentLi st .Target_Sequence . SegmentLi st . Launch 

Bu  rnout . LaunchAzDRDAl t . Down  rangeDi stance ' ) ; 

root . ExecuteCommand( 'Astrogator  */Satell i te/KVl  Setvalue 

Mai nsequence . SegmentLi st .Target_Sequence . SegmentLi st . Launch . Burnout . LaunchAzDRDAl t . Down rangeDi sta 
nee  B99.399  km'); 

root . ExecuteCommand( 'Astrogator  */Satel 1 i te/KVl  SetMCSControl Val ue 

Mai nsequence . SegmentLi st .Target_Sequence . Profi 1 es . Di fferenti al_Corrector  Launch 

Burnout. LaunchAzDRDAlt. DownrangeDi stance  Active  true'); 

root . ExecuteCommand( 'Astrogator  */Satel 1 i te/KVl  AddMCSSegmentControl 

Mai nsequence . SegmentLi st .Target_Sequence . SegmentLi st . Launch  Burnout . LaunchAzDRDAl t . LaunchAZ ' ) ; 
root . ExecuteCommand( 'Astrogator  */Satell i te/KVl  Setvalue 

Mai nsequence . SegmentLi  st .Target_Sequence . SegmentLi  st . Launch . Burnout . LaunchAzDRDAl  t . LaunchAz  - 
17.8115  deg'); 

root . ExecuteCommand( 'Astrogator  */Satel 1 i te/KVl  SetMCSControl Val ue 

Mai  nsequence . SegmentLi  st .Target_Sequence . Profi 1 es . Di fferenti al_Corrector  Launch 

Burnout. LaunchAzDRDAlt. LaunchAz  Active  true'); 

root . ExecuteCommand( 'Astrogator  */Satel 1 i te/KVl  AddMCSSegmentControl 

Mai nsequence . SegmentLi st .Target_Sequence . SegmentLi st . Launch  Launch . Epoch ' ) ; 

root . ExecuteCommand( 'Astrogator  */Satell i te/KVl  Setvalue 

MainSequence.SegmentList.Target_Sequence.SegmentList. Launch. Launch. Epoch  "1  Sep  2014  00:23:00" 
utcg ' ) ; 

root . ExecuteCommand( 'Astrogator  */Satel  1 i te/KVl  SetMCSControlVal ue 

Mai  nsequence .  SegmentLi  st .Target_Sequence . Profi 1 es . Di fferenti al_Corrector  Launch  Launch . Epoch 
Active  true') ; 

root . ExecuteCommand( 'Astrogator  */Satel 1 i te/KVl  AddMCSSegmentControl 

Mai nsequence. Segment Li st.Target_Sequence. Segment List. Maneuver  Impul si veMnvr. Cartesi an .X' ) ; 
root . ExecuteCommand( 'Astrogator  */Satel 1 i te/KVl  SetMCSControlVal ue 
Mai nsequence . SegmentLi st .Target_Sequence . Profi 1 es . Di fferenti al_Corrector  Maneuver 
impulsiveMnvr. Cartesian. X  Active  true'); 

root . ExecuteCommand( 'Astrogator  */Satel 1 i te/KVl  AddMCSSegmentControl 

Mai nsequence. Segment Li st.Target_Sequence. Segment List. Maneuver  Impul si veMnvr. Cartesi an. Y') ; 
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root . ExecuteCommand( 'Astrogator  */Satel 1 i te/KVl  SetMCSControl Val ue 

Mai nSequence . SegmentLi st .Target_Sequence . Profi 1 es . Di fferenti al_Corrector  Maneuver 

impulsiveMnvr. Cartesian. Y  Active  true'); 

root . ExecuteCommand( 'Astrogator  */Satei 1 i te/KVl  AddMCSSegmentControl 

Mai nSequence. Segment Li st. Target_Sequence. Segment Li st. Maneuver  ImpulsiveMnvr.Cartesian.Z') ; 
root . ExecuteCommand( 'Astrogator  */Satei 1 i te/KVl  SetMCSControl Val ue 
Mai nSequence . SegmentLi st .Target_Sequence . Profi 1 es . Di fferenti al_Corrector  Maneuver 
ImpulsiveMnvr.Cartesian.Z  Active  true'); 

root . ExecuteCommand( 'Astrogator  */Satel 1 i te/KVl  AddMCSSegmentControl 
Mai nSequence . SegmentLi st .Target_Sequence. SegmentLi  st . Propagate 
Stoppi ngcondi  ti ons . Du  rati  on . T ri pval ue ' ) ; 
root . ExecuteCommand( 'Astrogator  */Satell i te/KVl  Setvalue 

Mai nSequence . SegmentLi st .Target_Sequence. SegmentLi st . Propagate . Stoppi ngcondi tions . Du  rati  on .Tri  pva 
lue  135.336  sec'); 

root . ExecuteCommand( 'Astrogator  */Satel 1 i te/KVl  SetMCSControl Val ue 

Mai nSequence . SegmentLi st .Target_Sequence . Profi 1 es . Di fferenti al_Corrector  Propagate 

Stoppi  ngcondi  tions . Duration. TripValue  Active  true'); 


%  Configure  constraints 

root . ExecuteCommand( 'Astrogator  */Satell i te/KVl  Setvalue 

Mai nSequence . SegmentLi st .Target_Sequence. SegmentLi st . Propagate . Results  Epoch  X  Y  Z'); 
root . ExecuteCommand( [ 'Astrogator  */Satel 1 i te/KVl  SetMCSConstrai ntval ue 

Mai nSequence . SegmentLi st .Target_Sequence . Profi 1 es . Di fferenti al_Corrector  Propagate  X  Desired 

l 

»  ■  ■  ■ 

num2str(aim_poi nt(2)*le3)])  ; 

root . ExecuteCommand( 'Astrogator  */Satel  1  i te/KVl  SetMCSConstrai ntval ue 

Mai nSequence . SegmentLi st .Target_Sequence . Profi 1 es . Di fferenti al_Corrector  Propagate  X  Active 
true') ; 

root . ExecuteCommand( [ 'Astrogator  */Satel 1 i te/KVl  SetMCSConstrai ntval ue 

Mai  nSequence .  SegmentLi  st .Target_Sequence . Profi 1 es . Di fferenti al_Corrector  Propagate  Y  Desired 


num2str(aim_poi nt(3)*le3)]) ; 

root . ExecuteCommand( 'Astrogator  */Satel 1 i te/KVl  SetMCSConstrai ntval ue 

Mai nSequence . SegmentLi st .Target_Sequence . Profi 1 es . Di fferenti al_Corrector  Propagate  Y  Active 
true') ; 

root . ExecuteCommand( [ 'Astrogator  */Satel 1 i te/KVl  SetMCSConstrai ntval ue 

Mai  nSequence .  SegmentLi  st .Target_Sequence . Profi 1 es . Di fferenti al_Corrector  Propagate  Z  Desired 


num2str(aim_poi nt(4)*le3)]) ; 

root . ExecuteCommand( 'Astrogator  */Satel 1 i te/KVl  SetMCSConstrai ntval ue 

Mai nSequence . SegmentLi st .Target_Sequence . Profi 1 es . Di fferenti al_Corrector  Propagate  Z  Active 
true') ; 

[year, month, day, hour, minu, sec,  =  julian2greg(T_int) ; 

T_intUTC  =  [year  month  day  hour  minu  sec];  %  UTC  intercept  time 

T_intUTC  =  UTC2Str(T_i ntUTC) ; 

root . ExecuteCommand( [ 'Astrogator  */Satel 1 i te/KVl  SetMCSConstrai ntval ue 

Mai nSequence . SegmentLi st .Target_Sequence . Profi 1 es . Di fferenti al_Corrector  Propagate  Epoch  Desi red 


T_i ntUTC , '  UTCG ' ] ) ; 

root . ExecuteCommand( 'Astrogator  */Satel 1 i te/KVl  SetMCSConstrai ntval ue 

Mai  nSequence . SegmentLi st .Target_Sequence . Profi 1 es . Di fferenti al_Cor rector  Propagate  Epoch  Active 
true') ; 
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%  Initialize  corrector  sequence 

root . ExecuteCommand( 'Astrogator  */Satell i te/KVl  Setvalue 

Mai nsequence . SegmentLi st .Target_Sequence . Profi 1 es . Di fferenti al_Corrector . Mode  Iterate ' ) ; 
root . ExecuteCommand( 'Astrogator  */Satell i te/KVl  Setvalue 
MainSequence.SegmentList.Target_Sequence. Action  Run  active  profiles'); 

%  Run  MCS  sequence 

root . ExecuteCommand( 'Astrogator  */Satell i te/KVl  RunMCS'); 

root . ExecuteCommand( 'Astrogator  */Satel 1 i te/KVl  Appl yco r recti ons 

Mai nsequence . SegmentLi st .Target_Sequence ' ) ; 

root . ExecuteCommand( 'Astrogator  */Satel 1 i te/KVl  Cl earDWCGraphi cs ' ) ; 

root. ExecuteCommand(' Animate  *  Reset'); 

%  Extract  position  elements 

KVDP  =  KVl.DataProviders.ltem('Cartesian 

Position') . Group. Item( ' 12000 ' ) . Exec(scenario . starttime, scenario. stoptime, 1)  ; 

T  =  KVDP. DataSets. GetDataSetByName(' Time') .Getvalues; 

T  =  jul i andate(T) ; 

x  =  cell2mat(KVDP. DataSets. GetDataSetByName('x') .Getvalues) ; 
y  =  cell2mat(KVDP.DataSets.GetDataSetByName('y') .Getvalues) ; 
z  =  cell2mat(KVDP. DataSets. GetDataSetByName('z') .Getvalues) ; 

KVDP  =  KVl. DataProviders . Item( ' Cartesi  an 

Velocity') . Group. Item(' 12000') .Exec(scenario.starttime,scenario.stoptime,l) ; 
xdot  =  cell2mat(KVDP. DataSets. GetDataSetByName('x') .Getvalues) ; 
ydot  =  cell2mat(KVDP. DataSets. GetDataSetByName('y') .Getvalues) ; 
zdot  =  cell2mat(KVDP. DataSets. GetDataSetByName('z') .Getvalues) ; 

%  Save  Data 

KV. trajectory  =  [T, x,y,z, xdot, ydot, zdot] ; 


Determine  Intercept  Acquisition  Time 

This  step  calculates  the  time  at  which  the  KV  acquires  the  target  based  upon  the  global  intercept 
range  variable. 


r  =  inf;  %  initialize  range  to  target 

in  =  1;  %  loop  index 

while  r  >  int_aq 

R_KV  =  KV. trajectory(i n , 2 : 4) ; 

R  =  aim_point(2  :4)-R_KV; 
r  =  norm(R)  ; 
in  =  in+1; 

end 

KV.JD_int  =  KV.trajectory(in,l) ;  %  Julian  Day,  time  at  acquisition 

[year, month, day, hour, minu, sec, ~,~]  =  jul ian2greg(KV. JD_i  nt) ; 

KV.UTC_int  =  [year  month  day  hour  minu  sec]; 

end 
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Appendix  V:  TOM  Correlation  Algorithm  Script 


function  [KV,TOM_data]  =  TOMCorr(radar , KV, compl ex) 

%TOMCorr  processes  all  the  possible  combinations  of  TOMs  and  selects  the 

%arrangement  with  the  most  likely  probability. 

%  This  function  cannot  handle  mismatches.  If  a  plot_in  (plot  index)  is 
%  provided,  this  function  will  plot  all  the  possible  TOMs  on  a  subplot. 


Variable  Definition 

global  plotk; 

%  Subplot  Dimensions 
N  =  size(KV.TOM.l) ; 
a  =  cei 1 (N/2) ; 


TOM  Correlation 

Correlation  is  based  strictly  on  the  geometric  position  of  the  objects  in  the  focal  plane. 


%  %  Calculate  azimuth  and  elevation  from  radar  TOMs 
%  [~, radar_ang]  =  rangeangle(radar.TOM(: ,2:4) ') ;  %  rad 

%  %  output  in  form  [azimuth; elevation] 

%  radar_TOM  =  [(1:N)'  deg2rad(radar_ang(l, :)) '  deg2rad(radar_ang(2, :)) '] ; 
%  %  rad 

radar_TOM  =  radar. TOM2; 

%  All  possible  optical  TOM  arrangements 
optical_TOMs  =  cell(N,l); 

opti cal_TOM  =  kv.tom; 
opti  cal_TOMlDs  =  cell(N,l); 
opti cal _tomid  =  kv.tomid; 

TOM_actual  =  KV.TOM; 

r  =  zeros(N , 1) ; 

for  m  =  1 : N 

%  Save  previous  optical  TOM 
prev_optical  =  optical_TOM; 
prev_lD  =  opti cal _T0MID; 

%  Modify  previous  radar  TOM  to  next  arrangement 
if  m  >  1 

for  q  =  1 :  N 

in  =  q+1;  %  index  for  moving  TOM  entries 

if  in  >  si ze(prev_opti cal  ,  1) 
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%  correct  for  index  too  large 
in  =  i n-size(prev_optical  ,  1)  ; 

end 

opti cal_TOM(q , :)  =  prev_optical (i n , : ) ; 
opti cal_TOMlD{q , : }  =  prev_lD{i n , 1} ; 

end 

end 

%  Save  optical  TOM  arrangement 
opti cal_TOMs{m}  =  optical_TOM; 
opti cal_TOMlDs{m}  =  optical_TOMlD; 

%  Calculate  Correlation  for  particular  optical  TOM  arrrangement 
r(m)  =  corr2(radar_TOM( : , 2 : B) ,opti cal_TOMs{m} ( : , 2 : 3)) ; 

%  Plot  Each  Correlation 
if  ~isempty(plotk) 
figure(plotk) 
if  m  ==  1 
cl  f 

end 

subplot(a,2,m) ,  hold  on 

hi  =  plot(radar_TOM(: ,2) , radar_TOM(: , B) , 'vb') ; 
h2  =  plot(optical_TOM(:  ,2)  ,optical_TOM(:  ,3)  , 'Ar')  ; 
for  k  =  l:size(radar_TOM,l) 

plot([radar_TOM(k,2)  optical_TOM(k,2)] , [radar_TOM(k, 3) , . . . 
optical_TOM(k,3)] , 'Color' , [.2  k/size(radar_TOM,l)  0]); 

end 

title(['TOM  ' ,num2str(m) , ' ,  r=  ' ,num2str(r(m) ,4)]) 
xlabel ('Azimuth  (rad)') 
ylabel  ('Elevation  (rad)') 
axis  equal 

%1  egend([hl, h2] , ' Location ' , ' Eastoutsi de ' , ' Radar  Targets ' , . . . 

%'Optical  Targets') 

%  Label  Targets 

space  =  mean(mean([radar_TOM(: ,2)-optical_TOM(: ,2) , radar_TOM(: ,3)-. . . 

opti cal_TOM( :  ,3)]))  ;  %  Calculate  normalized  spacing 

for  p  =  1 : N 

text(radar_TOM(p, 2)+space , radar_TOM(p, 3)+space , num2str(radar_TOM(p, 1))) 
text (opti cal_TOM(p , 2) -space , opti cal_TOM(p , 3) -space , num2str (p)) 

end 

end 

end 

%  Create  additional  plot  space  if  needed 
%  i  f  m  >  a*2 
%  a  =  a+1; 

%  end 

%  subplot(a,2,m+l) ,  hold  on 

if  ~i sempty(pl otk) 
plotk  =  plotk+1; 

legend([hl  h2] , 'Location' , 'Best' , 'Radar  Targets ', 'Opti cal  Targets') 

end 
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%  Select  TOM  with  best  correlation 

%  r_max  is  the  max  correlation 
%  h  is  the  index  for  the  max  correlation 
[r_max,h]  =  max(r); 

TOM_data  =  [r_max,h,0];  %  Output  correlation  and  TOM  configuration 

%  Select  optical  TOM  of  particular  arrangment 
kv.tom  =  opti cal_TOMs{h} ; 
kv.tomid  =  opti cal_TOMiDs{h} ; 

%  Re-order  optical  charateri Stic  observations 

%  Original  data 

opti cal_char0  =  KV.char; 

mO  =  h-1;  %  shift  required 

for  k  =  1 : N 
m  =  mO; 
if  (k+m)  >  N 
m  =  m-N; 

end 

%  Update  optical  observations 
KV.char(k,:)  =  optical_charO(k+m, :) ; 

end 

%  Determine  if  correlation  was  successful 
space  =  max(radar_TOM( : ,2))/10; 
sue  =  1; 
for  k  =  1 : N 

%  if  ~strcmp(KV.TOMlD{k,l},complex{k,l}.lD) 

%  sue  =  0; 

%  end 

if  KV.TOM(k,l)  ~=  k 
sue  =  0; 

end 

end 

%  Plot  TOM 
if  ~i semptyCpl otk) 
figure(plotk) 
elf;  hold  on 
for  k  =  1 : N 

plot(radar_TOM(k,2) , radar_TOM(k, B) , 'Ar') 
plot(KV.TOM(k,2) ,KV.TOM(k,3) , 'vb') 
text(radar_TOM(k, 2)+space , radar_TOM(k , 3)+space,  . . . 
num2str(radar_TOM(k, 1))) 

text(TOM_actual (k,2)-space,TOM_actual (k,3)-space, . . . 

['A1 , num2str(TOM_actual  (k,l))]) 

plot( [radar_TOM(k, 2)  KV.TOM(k,2)] , [radar_TOM(k, 3)  KV.TOM(k,3)] 
'Color', [.2  k/size(radar_TOM,l)  0]); 
xlabel ('Azimuth  (rad)') 
ylabel  ('Elevation  (rad)') 

title(['Selected  TOM  Configuration,  r  =  ' , num2str(r_max, 3)] ) 
plotk  =  plotk+1; 

end 

end 
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T0M_data(3)  =  sue; 
end 
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Appendix  VI:  Dempster-Shafer  Discrimination  Algorithm  Coding 


The  code  for  the  Dempster-Shafer  algorithm  was  developed  by  Dr.  Sami  Amashi,  MIT- 
Lincoln  Laboratory,  Group  3 1 . 

function  [PPs]  =  DempsterShaferOpen(means, sigmas, measures) 

%DempsterShaferOpen  calculates  the  pignistic  probabilities  of  a  given 
%series  of  M  measurements  belonging  to  each  class  given  a  set  of  observations 
%and  the  expected  distributions  of  those  measurements  for  each  class.  It 
%utilizes  an  open-world  assumption  and  does  not  eliminate  any  conflict  in 
%the  measurements.  N  is  the  number  of  objects  observed 

%  Inputs:  means  -  expected  observations,  expected  as  a  MxN  row  vector 

%  sigmas  -  observation  standard  deviations  as  a  lxN  row  vector 

%  obv  -  measurements  expected  as  an  Mx(N+l)  matrix,  where 

%  the  first  row  is  object  indices. 

% 

%  Outputs:  PPs  -  Pignistic  Probabilities  for  each  class. 

% 

%  The  Dempster-Shafer  functionality  of  this  program  was  writen  by  Sami 
%  Amashi,  MIT-Lincoln  Laboratories.  This  function  merely  assembles  his 
%  functions  into  one  function. 

global  num_class; 

N  =  size(measures , 1) ;  %  Number  of  objects 

M  =  size(means , 1) ;  %  Number  of  characteristics 

PPs  =  zeros(N , num_cl ass) ; 

for  k  =  1:N  %  Loop  over  each  detection 

for  m  =  2 : (M+l)  %  Loop  over  each  charateri sti c 

obv  =  measures(k,m) ;  %  observed  characteristic 

%  Compute  probablity  of  each  measurement  from  each  distribution 

mus  =  means(m-l, :) ;  %  expected  observations 

sigma  =  sigmas(m-l);  %  observation  SD 

[pbl, pb2 , pbB]  =  ProbAssign3(obv,mus, sigma) ; 

class_prob  =  [pbl;  pb2 ;  pbB]  ;  %  probability  of  each  observation  from  each  class 

%  Assign  evidence  measure 

EM  =  DS_Consonant_UnNormalized(class_prob) ; 

%  Combine  previous  and  current  evidence  measures 
if  m  >  2  %  skip  combination  for  first  characteristic 
EM_prev  =  DS_OWCombi nation (EM, EM_prev) ; 

el  se 

EM_prev  =  EM; 

end 

end 

%  Calculate  pignistic  probabilities  based  on  open  world 
PPs(k,:)  =  DS_PignisticProbability(EM_prev,l) ' ; 

end 
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end 
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